text:
# Algorithms of maximum likelihood data clustering with applications

## 1 Introduction

One aspect of the Information Technology revolution is that huge amounts of data have become available. For example, every transaction in many financial markets is recorded and bio-informatics technology allows us today to monitor genome wide gene expression. These informations may allow for a quite detailed description of such complex systems as a financial market or a cell, to quote just two examples. However our understanding of a complex systems is, in many cases, limited by our ability to efficiently organize massive streams of information. As a result, methods for handling, organizing or mining large data sets have become of great practical importance.

Data clustering deals with the problem of classifying a set of $`N`$ objects into groups so that objects within the same group are more similar than objects belonging to different groups. Each object is identified by a number $`D`$ of measurable features: Hence object $`i=1,\mathrm{},N`$ can be represented as a point $`\stackrel{}{x}_i=(x_i^{(1)},\mathrm{},x_i^{(D)})`$ in a $`D`$ dimensional space. Data clustering aims at identifying clusters as more densely populated regions in this vector space. More precisely, a configuration of clusters is represented by a set $`𝒮=\{s_1,\mathrm{},s_N\}`$ of integer labels, where $`s_i`$ is the cluster to which object $`i`$ belongs.

We focus on unsupervised approaches, which do not require any further information apart from that contained in the data. These are of great practical importance specially in the presence of huge data sets. We assume, in other words, that the data set is homogeneous: all features are equally important or are already appropriately weighted.

The classic approaches to data clustering are partitioning methods and hierarchical clustering. Partitioning methods are based on two elements: 1) a distance between objects, which allows one to measure their similarity and 2) a cost function whose minima correspond to “optimal” clustering configurations. For example, a typical K-means (KM) approach takes as cost function the sum of squared distances of objects to the centroid $`\stackrel{}{X}_s`$ of the cluster in which they are classified :

$$H_{KM}\{𝒮\}=\underset{s}{}\underset{i:s_i=s}{}\stackrel{}{x}_i\stackrel{}{X}_s^2,\stackrel{}{X}_s=\frac{1}{n_s}\underset{i:s_i=s}{}\stackrel{}{x}_i$$
(1)
where $`n_s=_{i:s_i=s}1`$ is the number of objects in cluster $`s`$, and $`\stackrel{}{x}^2=_j(x^{(j)})^2`$ defines the Euclidean distance.

Apart from cases where a cost function is naturally suggested by the problem itself111The typical example is that of deciding the location of $`K`$ distribution centers which should serve $`N`$ cities. It is natural to look for the locations which minimize the sum of distances., the choice of the cost function or of a distance is an element of arbitrariness. A further problem of these approaches is that one needs to predefine the number $`K`$ of clusters from the beginning. Note, for example that $`H_{KM}`$ of Eq. (1) attains its minimal value $`H_{KM}=0`$ when each object is in a different cluster ($`n_s=1`$). The number $`K`$ of clusters should then be fixed a priori or one has to introduce a $`K`$-dependent term in the cost function. This is a further element of arbitrariness.

The problem of finding the cluster structure which minimizes the cost function may be quite difficult: depending on the form of the cost function and on the data set the ”cost landscape” can be either simple, with a single minimum which is easily accessible dynamically, or complex, with many metastable minima. Data clustering methods also differ for the specific algorithms used to reach a (local) minimum.

A second, very popular approach to data clustering – called hierarchical clustering – is based on the definition of a distance between objects and clusters of objects and a very simple algorithm: Given a configuration with $`K>1`$ clusters, it merges the two closest clusters into a single one. In this way, starting from the configuration with $`K=N`$ clusters, the algorithm generates a sequence of configurations as $`K`$ varies from $`N`$ to $`1`$. This sequence of configurations and their hierarchic organization, can be represented by a convenient and compact graphical tool called dendrogram . The minimal spanning tree algorithm – also called Single Linkage (SL) – the Average Linkage (AL) and the Centroid Linkage (CL) algorithms are examples of hierarchical clustering methods222These algorithms differ in the way the distance between clusters is computed: In SL (AL) the distance between clusters $`s`$ and $`r`$ is the minimal (average) distance between items in cluster $`s`$ and items in cluster $`r`$; in CL the distance $`\stackrel{}{X}_s\stackrel{}{X}_r`$ between centroids is used.. Rather than a single cluster structure, this algorithm provides a hierarchic sequence of cluster structures. The choice of the best cluster structure is left arbitrary. Applications of this approach to real world data is discussed in Refs. .

Many other approaches to data clustering have been proposed: For example Refs.  used singular value decomposition to identify clusters. Identifying principal components with cluster structures imposes an orthogonality constraint between clusters which may be unnatural. Expectation Minimization is a further approach  where the density of points is modeled as a mixture of Gaussians whose centering and scale parameters are fit by maximum likelihood. Kohonen et al.  have proposed an algorithm based on Self-Organizing Maps (SOM) whereas Blatt et al. proposed the Super-Paramagnetic Clustering (SPC) method . The latter is based on a mapping to an interacting particle system whose “magnetic” properties describe the cluster structure of data and has been applied to a range of problems (see e.g. ). These methods rely on ad hoc definitions of a dynamics (SOM) or of the particle-particle interaction (SPC) which are tuned by several parameters or functions.

In the following sections we briefly review yet a different approach to dataclustering that we recently proposed , we describe possible algorithms based on it, and we finally compare the outcomes with those of standard data clustering methods, as those described above.

## 2 Maximum likelihood data clustering

We have devised  a fully unsupervised, parameter free approach to data clustering which derives from a maximum likelihood (ML) principle. The key idea is that objects are similar if they have something in common. In a correct classification, objects belonging to the same cluster should share a common component:

$$\stackrel{}{x}_i=g_{s_i}\stackrel{}{\eta }_{s_i}+\sqrt{1g_{s_i}^2}\stackrel{}{ϵ}_i.$$
(2)
Here $`\stackrel{}{x}_i`$ is the vector of features of object $`i`$, normalized so that $`_tx_i^{(t)}=0`$ and $`\stackrel{}{x}_i^2=_t[x_i^{(t)}]^2=D`$ for all $`i=1,\mathrm{},N`$333Equivalently, one may generalize Eq. (2) by adding a constant term $`r_i\stackrel{}{1}`$ and a scale factor $`\sigma _i`$. The maximum likelihood estimates of these parameters are the mean and the variance. Subtracting the mean and rescaling by the variance, leaves us with the normalized data set. The parameters $`r_i`$ and $`\sigma _i`$ are irrelevant as far as the cluster structure is concerned. The latter indeed only depends on the “internal structure” of correlations.; $`s_i`$ is the label of the cluster to which it belongs. $`\stackrel{}{\eta }_s`$ is the vector of features of cluster $`s`$ and $`g_s`$ tunes the similarity of objects within cluster $`s`$: For $`g_s=1`$ all objects with $`s_i=s`$ are identical whereas when $`g_s`$ is small objects are very different. The cluster index $`s`$ ranges from $`1`$ to $`N`$ in order to allow also for the case of $`N`$ clusters of one object each. $`\stackrel{}{ϵ}_i`$ describes the deviation of the features of object $`i`$ from the cluster’s features and measurement errors. We take Eq. (2) as a statistical hypothesis and assume that both $`\stackrel{}{\eta }_s`$ and $`\stackrel{}{ϵ}_i`$ (for all $`i,s=1,\mathrm{},N`$) are Gaussian vectors with zero mean and variance $`E[(\eta _s^{(t)})^2]=E[(ϵ_i^{(t)})^2]=1`$. For any given set of parameters $`(𝒢,𝒮)=(\{g_s\},\{s_i\})`$ it is possible  to compute the probability (density) $`P(\{\stackrel{}{x}_i\}|𝒢,𝒮)`$ of observing the data set $`\{\stackrel{}{x}_i\}`$ as a realization of Eq. (2). From this it is possible  to compute the likelihood444Note that an assumption on the prior probability $`P(𝒢,𝒮)`$ is invoked by Bayes formula in this passage. We take $`P(𝒢,𝒮)=\text{const}`$ which means that every cluster structure $`𝒮=\{s_i\}`$ is a priori equiprobable. $`P(𝒢,𝒮|\{\stackrel{}{x}_i\})`$. It turns out  that the resulting expression depends only on the Pearson’s coefficient

$$C_{i,j}=\frac{\stackrel{}{x}_i\stackrel{}{x}_j}{\sqrt{\stackrel{}{x}_i^2\stackrel{}{x}_j^2}}.$$
(3)
For a given cluster structure $`𝒮`$, the likelihood is maximal when the parameters $`g_s`$ take the values

$$g_s^{}=\sqrt{\frac{c_sn_s}{n_s^2n_s}}$$
(4)
if $`n_s>1`$ and $`g_s^{}=0`$ if $`n_s1`$. Here $`n_s`$ is the number of objects in cluster $`s`$ and

$$c_s=\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}C_{i,j}\delta _{s_i,s}\delta _{s_j,s}.$$
(5)
Note indeed that when $`c_sn_s`$, as for uncorrelated objects, then $`g_s^{}0`$ whereas if objects are very similar, $`c_sn_s^2`$ and $`g_s^{}1`$.

The maximum likelihood of structure $`𝒮`$ can be written as $`P(𝒢^{},𝒮|\{\stackrel{}{x}_i\})e^{D_c(𝒮)}`$, where the log-likelihood per feature $`_c`$ is given by

$$_c(𝒮)=\frac{1}{2}\underset{s:n_s>1}{}\left[\mathrm{log}\frac{n_s}{c_s}+(n_s1)\mathrm{log}\frac{n_s^2n_s}{n_s^2c_s}\right].$$
(6)
This depends on the original data through the coefficients $`c_s`$ of Eq. (5). The function $`_c`$ provides a likelihood measure for cluster structures. The ML structure is that which maximizes $`_c`$.

There are several interesting features of $`_c`$:

* If the objects are unrelated ($`g_s^{}=0`$ or $`c_s=n_s`$) or if they are classified in singleton clusters ($`n_s=1`$ for all $`s`$) we find $`_c=0`$. Loosely speaking, $`\mathrm{max}_𝒮_c(𝒮)`$ measures the amount of structure present in the data-set.
* The maxima of $`_c`$ do not necessarily coincide with a single cluster containing all objects – as in the SPC approach of Ref.  – nor with the configuration with all objects in different clusters – as for $`H_{KM}`$555Note that Eq. (1) can be written as $`H_{KM}=2_{s:n_s>0}(n_sc_s/n_s)`$.
* $`_c`$ does not depend on any parameter.
* The number $`K`$ of clusters is not fixed a priori. Rather it is predicted.
* The interpretation of the results is transparent in terms of the model (2).

Eq. (2) may not be the most appropriate description for a particular data set. In much the same way, a straight line may not be the best description of a set of points on a plane. Still in this example, least squares provide an unambiguous method to compute the coefficients and a statistical measure of the goodness of fit. The same is true for our method: A sharp maximum of $`_c/N`$ indicates a robust cluster structure whereas when $`_c/N`$ is small and has several local maxima the cluster structure is not statistically significant.

The statistical hypothesis Eq. (2) is specially helpful for high dimensional data sets ($`D1`$) where geometric intuition becomes problematic. From the point of view of the computational cost, the dimensionality $`D`$ of the data set enters only in the calculation of the matrix $`C_{i,j}`$. For small $`D`$ methods based on geometry and visual inspection may be preferrable.

Our approach has some similarity with Expectation Maximization data clustering , which is also based on likelihood maximization. However the statistical hypothesis in Ref.  is very different from Eq. (2).

## 3 Data clustering algorithms

The function $`_c(𝒮)`$ is a measure of likelihood for cluster structures and it can be used to compare cluster structures produced by standard methods. For example, it allows one to compare classifications with a different number of clusters and hence to select the optimal number of clusters in KM, SL, AL or CL algorithms. We shall compare later the cluster structure found in this way with those found with algorithms based on $`_c`$ itself. Most importantly in fact, $`_c`$ can be used as the basis for clustering algorithms.

A very powerful method to find ML configurations is simulated annealing (SA)  with cost function $`_c`$. This is simply implemented via Metropolis dynamics  on the dynamical variables $`s_i`$ with progressively decreasing (fictitious) temperature $`T`$. This method produces a cluster configuration, which, if the annealing schedule is appropriately chosen, has a good chance of being that of maximum likelihood. The behavior of the system as a function of the fictitious temperature is similar to that described in Ref. . However, the SPC method (in its original formulation ) yields a trivial configuration with a single cluster in the limit $`T0`$. This is a consequence of the mapping of the data set into a particle system with ferromagnetic interactions, corresponding to positive correlations. This forces one to study intermediate temperatures666Modifications to take into account negative correlations are suggested in Ref. .. The ML configuration found in our method when $`T0`$ is already non-trivial. Further insight on the cluster structure can be obtained by studying the system at finite $`T`$, as shown in Ref. . For simplicity however, we will limit the present discussion to the $`T0`$ limit.

A second possibility is that of using a deterministic maximization (DM) technique. We discuss below a very simple algorithm: Given a configuration $`S`$, for all $`i=1,\mathrm{},N`$ (in some fixed order) propose all moves $`s_is`$ and compute the corresponding variation $`\mathrm{\Delta }_{s_is}`$ of $`_c`$; Find $`s_i^{}=\mathrm{arg}\mathrm{max}_s\mathrm{\Delta }_{s_is}`$ and set $`s_i=s_i^{}`$. Repeat these steps until no change occurs. This will lead in the end to a local maximum of $`_c`$ which we shall call the $`𝒮_{DM}`$ configuration.

Finally we propose a merging algorithm (MR) in the spirit of hierarchical clustering based on $`_c`$.

1. Start from $`N`$ clusters composed of one object each (e.g. $`s_i=i`$ $`i=1,\mathrm{},N`$).
2. At each step of the algorithm merge two clusters into a single one in such a way that the cost $`_c`$ of the resulting configuration is minimal.
3. Repeat step 1) $`N1`$ times until the configuration with a single cluster is reached.

Let $`𝒮_{MR}`$ be the state with maximal $`_c`$ found with this algorithm. $`𝒮_{MR}`$ is not necessarily a local maximum of $`_c`$ with respect to single “spin-flip” moves $`s_is_i^{}`$, while $`𝒮_{SA}`$ and $`𝒮_{DM}`$ are local maxima. Hence, we expect that $`_c(𝒮_{MR})_c(𝒮_{DM})_c(𝒮_{SA})`$. On the other hand, the hierarchical clustering algorithm is definitely faster than both of the others: Even if it delivers only an approximation to the maximum likelihood, faster execution can be crucial when dealing with huge data sets.

In addition, the hierarchical clustering technique MR offers a very convenient graphical representation of the cluster structure in terms of dendrograms. Usually dendrograms report the tree structure as a function of a resolution length. As this length increases, clusters are merged one by one until a single cluster remains when the resolution exceeds the maximal distance $`\mathrm{max}_{i,j}\stackrel{}{x}_i\stackrel{}{x}_j`$ between data points. We draw the dendrogram as a function of the log-likelihood per feature $`_c`$. At each merging step $`s,rq`$ let $`\mathrm{}_s`$, $`\mathrm{}_r`$ be the contribution of clusters $`s`$ and $`r`$ respectively to $`_c`$ before merging and $`\mathrm{}_q`$ be the corresponding cluster log-likelihood of the merged cluster $`q`$. In the dendrogram, we merge two points at height $`\mathrm{}_s`$ and $`\mathrm{}_r`$ into a single point at height $`\mathrm{}_q`$. There are three types of branching points: i) if $`\mathrm{}_q>\mathrm{}_s+\mathrm{}_r`$ the step increases the likelihood of the configuration ii) if $`\mathrm{}_q<\mathrm{}_s+\mathrm{}_r`$ but $`\mathrm{}_q>\mathrm{max}(\mathrm{}_s,\mathrm{}_r)`$ the global likelihood decreases but the cluster likelihood increases. Finally iii) steps with $`\mathrm{}_q<\mathrm{max}(\mathrm{}_s,\mathrm{}_r)`$ represent unlikely merging moves.

## 4 Applications and discussion

We discuss the application of the above ideas to two different data sets: The first consists of daily returns of the NYSE from January the 1st 1987 to March the 30th 1999. We consider $`N=1000`$ of the assets traded most frequently throughout this period . If $`p_i(t)`$ is the opening price of asset $`i`$ in day $`t`$, we set

$$x_i^{(t)}=\frac{\mathrm{log}p_i(t)/p_i(t1)r_i}{\sigma _i}$$
(7)
where $`r_i`$ and $`\sigma _i`$ enforce normalization ($`\stackrel{}{x}_i\stackrel{}{1}=0`$ and $`\stackrel{}{x}_i^2=D`$). Actually it is convenient to introduce a further linear transformation in order to eliminate the so called “market mode”: $`\stackrel{}{x}_i^{}=\alpha _i(\stackrel{}{x}_i\stackrel{}{x}_{\mathrm{avg}})`$ where $`\stackrel{}{x}_{\mathrm{avg}}=_i\stackrel{}{x}_i/N`$ so that $`_i\stackrel{}{x}_i^{}=0`$ ($`\alpha _i`$ is a coefficient to restore $`\stackrel{}{x}_i^{}^2=1`$).

Secondly we discuss the gene expression data set measured in  and already studied in Refs.  with a variety of techniques. The data measure the expression of the 2467 genes of known function of the yeast Saccharomyces cerevisiæ, whose complete genome has been sequenced. In particular, we focused on the 18 time values coming from the alpha-factor block-and-release experiment in ; these correspond to approximately two cell cycles. Also in this case, we take the logarithm of the value of the gene expression, so that the data distribution is closer to a Gaussian, and normalize the sets as required by Eq. 2.

In both cases the data set is homogeneous by definition: features characterize the system – the market or the cell – at different times and there is no a priori reason to weight features differently777Note however that some major crashes such as black Monday of October 1987 are included in the financial data set.. The main objective of data clustering for these data sets is that of identifying groups of assets with a similar market behavior or groups of genes with a similar function. This information is of considerable interest for risk management strategies and for the understanding of biological functions respectively in the two cases. We shall not enter here into details of either finance or genetics, but rather focus on the properties of data clustering algorithms for these two case studies. On one hand we are interested in identifying fast clustering algorithms, on the other it is desirable that the resulting cluster structures does not depend strongly on the detailed algorithm used. In few words, one would like to identify fast algorithms which provide results which are consistent with more elaborate and time consuming ones. Note that the two data sets are very different in nature: For market data time series are very long ($`D=3114`$) and correlations are very weak while gene expression is recorded for few ($`D=18`$) time points but correlations are very strong. Thus they constitute a good testing ground for our analysis.

The cluster structure which emerges from the minimization of $`_c`$ are quite meaningful: For market data we find a strong overlap with the classification of stocks in sectors of economic activity. In the gene expression data set, as in Ref. , our method identifies clusters related to biological functions. This information is preserved more or less in the ML cluster structures found with different methods.

Standard hierarchical clustering methods (SL, CL and AL) give ML configurations which differ considerably from that found with other methods (see Ref.  for details). Table 1 shows that the configurations $`𝒮_{SL}`$, $`𝒮_{AL}`$ and $`𝒮_{CL}`$ have many more clusters than the others and a much lower value of $`_c`$. Fig. 1 compares the log-likelihood of configurations produced in different algorithms as a function of the number of clusters. The curve corresponding to AL deviates markedly from those of other methods; the SL and CL curve (not shown) deviate even more. Not only ML configurations are very different because of a larger number of clusters. Even configurations with a fixed number $`K`$ of clusters produced by the SL, CL and AL algorithm differ considerably among themselves and from the corresponding configurations with $`K`$ clusters found with other methods (see Fig. 2).

The closeness of points MR, DM and SA in Fig. 1 suggests that methods based on $`_c`$ produce consistent results, irrespective of the algorithm. The closeness of two cluster structures $`𝒮`$ and $`𝒮^{}`$ can be quantified in terms of overlaps: Consider the intersection $`𝒮𝒮^{}`$ of the two cluster structures – defined so that objects $`i`$ and $`j`$ are in the same cluster in $`𝒮𝒮^{}`$ only if they are in the the same cluster both in $`𝒮`$ and in $`𝒮^{}`$. We define a geometric overlap

$$O_g(𝒮|𝒮^{})=\frac{|𝒮𝒮^{}|}{|𝒮^{}|}$$

where $`|𝒮|`$ is the number of pairs of objects that are in the same cluster in configuration $`𝒮`$. In words, $`O_g(𝒮|𝒮^{})`$ is the fraction of pairs of objects that are in the same cluster in $`𝒮^{}`$ which also belong to the same cluster in $`𝒮`$. Equivalently, it is the probability that that two randomly chosen objects which belong to the same cluster in $`𝒮^{}`$ are also in the same cluster in $`𝒮`$. Tables 2 and 3 report the geometric overlaps between the ML configurations and confirm the conclusion that methods based on $`_c`$ yield consistent results. In fact, the overlaps among configurations MR, DM and SA are generally larger than those in the sector SL, AL, KM. These features manifest strongly in weakly correlated data sets, such as market data, and to a milder extent when correlations are stronger, as in the gene expression data set. Still the distribution of cluster sizes produced by SL and AL algorithms is much more uneven than that found with other methods.

The geometric overlap is dominated by clusters which contribute negligibly to the likelihood and which are therefore not statistically significant. It is then preferable to introduce a likelihood overlap

$$O_l(𝒮^{}|𝒮)=\frac{_c(𝒮𝒮^{})}{_c(𝒮)}$$

which yields the fraction of likelihood of configuration $`𝒮`$ which is also “explained” by configuration $`𝒮^{}`$. Tables 4 and 5 show that most of the relevant information (in terms of likelihood) gained by standard methods is also contained in ML configurations found with our method. The opposite is however not true.

The discrepancy between results of deterministic minimization and hierarchical clustering algorithms is much reduced if one uses $`_c`$. Indeed, MR and DM predict a number of clusters $`K`$ of the same order whereas AL predicts a much larger $`K`$ than KM, specially if correlations are weak, as in the financial data set (see table 1). As a further evidence, Fig. 2 compares the overlap between configurations AL and KM with $`K`$ clusters with the overlap between the configuration MR with $`K`$ clusters and the corresponding configuration DM (see caption): Configurations MR and DM are more similar than configurations AL and KM.

At odds with SL, AL and CL algorithms, hierarchical clustering yields a ML configurations (MR) which is quite consistent with maximization algorithms (SA or DM) which are computationally more demanding. This, we believe, is due to the fact that $`_c`$ generates in typical cases a configuration landscape which is smoother and with more easily accessible maxima, than that generated by conventional distance-based cost functions, such as $`H_{KM}(𝒮)`$. This conjecture is partially confirmed by looking at simple cases. Take, for example, a data set generated by Eq. (2) with constant correlation $`\stackrel{}{x}_i\stackrel{}{x}_j=g>0`$ for all $`ij`$. All configurations with two clusters have the same cost $`H_{KM}`$, and there are $`2^{N1}1`$ such configurations. But the configurations of maximal $`_c`$ with two clusters are only those with an isolated object and a large cluster of size $`N1`$ (there are $`N`$ such configurations). Thus $`_c`$ gives a much smaller degeneracy than $`H_{KM}`$ in this case.

The MR algorithm delivers information on the hierarchical organization of the data set. The dendrograms for the two case studies are shown in Figs. 3 and 4. The dendrogram of market data reveals a complex structure consistent with the scaling laws reported in Ref.  (see inset). A much more uniform cluster distribution is found in gene expression data. Figs. 3 and 4 only report the first two types of branch points: i) those which increase the total likelihood and ii) those which do not increase the total likelihood but where the likelihood of the unified cluster is larger than those of both individual clusters. These two sets of branching points are clearly separate and describe the cluster structure at different resolutions. In most branching of type i) a large cluster merges with a single object whereas type ii) steps merge clusters into larger clusters.

Gene expression data shows a rich hierarchical structure of type ii) branchings. The ML structure has a large number of clusters with $`8÷15`$ elements. This suggests that the model Eq. (2) does not provide an exhaustive description of the correlations of this data set. There are indeed sizeable correlations between clusters that are not taken into account. An alternative way to describe these correlations, besides that provided by the MR algorithm, is to perform a reclustering of the cluster patterns $`\stackrel{}{X}_s`$ with the SA algorithm. More precisely this amounts to taking $`\stackrel{}{x}_s=\stackrel{}{X}_s/\stackrel{}{X}_s`$ as the input data set. Reclustering steps can be iterated until a configuration with singleton clusters is found. This leaves us with a sequence of nested reclustered configurations which provide a hierarchical description of the correlations present in the original data set. In fig. 5 we show the effect of a second application of the SA clustering procedure on the biological data set.

## 5 Conclusions

We have discussed a maximum likelihood approach to data clustering and its application to two examples. The method is based on a simple statistical description of data where similar objects have something in common (Eq. 2). The likelihood that a particular data set is described by such a model with a given cluster structure can be efficiently computed. This provides a parameter-free measure for cluster structures which can then be taken as the basis of clustering algorithms.

Different algorithms to find maxima of the likelihood – or good approximations to them – have been introduced. On one hand we discuss computationally expensive methods, such as simulated annealing technique (SA) and deterministic minimization (DM) algorithms, which are expected to provide a good approximation to the maximum likelihood structure. At odd with most other methods, the number of clusters should not be fixed a priori but is predicted by our algorithms. On the other hand, we introduce a simple deterministic hierarchical clustering algorithm (MR) which provides a much faster approximation.

We show that standard methods with distance based measures of clustering predict cluster structures which may be very different according to the specific algorithm used (SL, AL or KM). On the contrary, the maximum likelihood structures predicted by different algorithms based on the measure $`_c`$ (MR, DM and SA) are very similar. This suggests that, for instances where a meaningful cluster structure exists, the log-likelihood landscape has a broad and easily accessible maximum, whereas distance-based cost functions such as $`H_{KM}`$ (Eq. 1) may give rise to a more complex landscape.

Our approach is particularly suited for large dimensional data sets. Indeed, the dimensionality $`D`$ only enters in the calculation of Pearson’s coefficients $`C_{i,j}`$.

Data clustering has been regarded as an ill defined problem (see e.g. ). Indeed in conventional approaches one needs to define a similarity measure and a resolution scale. The approach proposed in the present work does not suffer from these drawbacks. The cluster structure is entirely determined by the internal correlations of the data set. This makes data clustering a well defined problem.

We are grateful to R. N. Mantegna for providing the financial data set.


format: text

text:
# A New Proof of the Faithfulness of Brauer’s Representation of Temperley-Lieb Algebras

## 1 Introduction

Temperley-Lieb algebras, whose roots are in statistical mechanics (see  and , Appendix II.b), play a prominent role in knot theory and low-dimensional topology. They have entered this field through Jones’ representation of Artin’s braid groups in a version of these algebras (see , §11, and works of Jones cited therein), on which the famous Jones polynomial of knot theory is based (see also , , ). The version of Temperley-Lieb algebra we deal with here, which is pretty standard (see , and works by Kauffman and others cited therein), is defined as follows.

The Temperley-Lieb algebra $`𝒯_n`$, where $`n2`$, is an associative algebra over a field of scalars whose nature is not important for us in this paper (in , p. 9, it is taken to be the field of rational functions $`P/Q`$ with $`P,Q\text{Z}[x,x^1]`$). The basis of $`𝒯_n`$, which we will call $`𝒥_n`$, is freely generated from a set of generators $`\text{1}^n,h_1^n,\mathrm{},h_{n1}^n`$ with multiplication subject to the monoid equations (i.e., multiplication is associative and $`\text{1}^n`$ is a unit for it) and the equations

$$\begin{array}{cc}\text{(1)}\hfill & h_i^nh_j^n=h_j^nh_i^n\text{}ij2\text{,}\hfill \\ (2)\hfill & h_i^nh_{i\pm 1}^nh_i^n=h_i^n,\hfill \\ (3)\hfill & h_i^nh_i^n=h_i^n.\hfill \end{array}$$

For the definition of multiplication in $`𝒯_n`$ the equation $`(3)`$ is replaced by

$$\begin{array}{cc}(3p)\hfill & h_i^nh_i^n=ph_i^n\hfill \end{array}$$

with $`p`$ a specified nonzero scalar (in , it is taken to be $`x^2x^2`$, but other values are found in other treatments of Temperley-Lieb algebras). The basis $`𝒥_n`$ of $`𝒯_n`$ is finite: its cardinality is the $`n`$-th Catalan number $`(2n)!/(n!(n+1)!)`$. The algebra $`𝒯_n`$ is the vector space whose basis is $`𝒥_n`$.

For those $`𝒯_n`$ where $`p`$ in $`(3p)`$ is a natural number greater than or equal to 2 there is a representation in matrices due to Brauer (see  and , Section 3). This is the orthogonal group case of Brauer’s representation restricted to the Temperley-Lieb subalgebra of the Brauer algebra. The faithfulness, i.e. isomorphism, of this representation was established by Jones in  (Section 3, Theorem 3.4, p. 330) by referring to his technique of Markov trace. Our purpose in this paper is to give a different, elementary and self-contained, proof of that faithfulness. We believe our proof is worth publishing because of its aesthetic value.

The computational interest of Brauer’s representation is lessened by the fact that for $`𝒯_n`$ we have to pass to $`p^n\times p^n`$ matrices, and so get an exponential growth. However, this growth is to be expected, since for sufficiently large $`n`$ the $`n`$-th Catalan number is greater than $`(2\epsilon )^{2n}`$ (this can be computed with the help of Stirling’s Theorem). When they were first introduced in , Temperley-Lieb algebras were represented in matrices in a manner different from Brauer’s (see , Appendix II.b, p. 264), again with an exponentional growth.

Brauer’s representation is provided by assigning to the elements of the basis $`𝒥_n`$ of $`𝒯_n`$ particular 0-1 matrices that satisfy (1), (2) and (3$`p`$) for $`p`$ a natural number greater than or equal 2. For the faithfulness of the representation it suffices to establish that the list of matrices assigned to the elements of $`𝒥_n`$ is linearly independent. Our proof of linear independence proceeds as follows. We introduce a linear order on the matrices, and establish that every matrix has an entry with 1 where all the matrices preceding it in the order have 0.

It is shown in  that Brauer’s representation of Temperley-Lieb algebras is based on the fact that the Kronecker product of matrices gives rise to an endofunctor of the category of matrices that is adjoint to itself. The result of this paper is interesting because of applications in the area covered by .

## 2 Ordering $`𝒥_n`$

For $`n2`$ and $`1in1`$, let $`h_{i,1}^n`$ be $`h_i^n`$, and for $`2ki`$ let $`h_{i,k}^n`$ be $`h_{i,k1}^nh_{ik+1}^n`$. We call the expressions $`h_{i,k}^n`$, for $`1ki`$, blocks.

Every element of $`𝒥_n`$ different from $`\text{1}^n`$ is denoted by a unique expression in Jones’ normal form

$$h_{i_1,k_1}^n\mathrm{}h_{i_l,k_l}^n$$

where $`l1`$, $`i_1<\mathrm{}<i_l`$ and $`i_1k_1<\mathrm{}<i_lk_l`$ (see , §4.1.4, , §2.8,  and , §10). Let $`h_{i,k}^n`$ stand for $`\text{1}^n`$ when $`k=0`$, and let us call these expressions too blocks. Then every element of $`𝒥_n`$ is denoted by a unique expression in the normal form

$$h_{1,k_1}^n\mathrm{}h_{n1,k_{n1}}^n$$

where $`0k_ii`$, and for every $`k_i,k_j>0`$ such that $`i<j`$ we have $`ik_i<jk_j`$. So we may identify every element of $`𝒥_n`$ by a sequence $`k_1\mathrm{}k_{n1}`$ subject to the conditions on $`k_i`$ and $`k_j`$ we have just stated.

Let these sequences be ordered lexicographically from the right; in other words, let

$$k_1^{}\mathrm{}k_{n1}^{}<k_1\mathrm{}k_{n1}$$

if and only if for some $`i`$ we have $`k_i^{}<k_i`$ and for every $`j>i`$ we have $`k_j^{}=k_j`$. This induces a linear order on $`𝒥_n`$, which we call the lexicographical order of $`𝒥_n`$.

## 3 Matrices, relations and graphs for $`𝒥_n`$

Let $`p`$ be a natural number greater than or equal to 2, and let $`E_p`$ be the $`1\times p^2`$ matrix that for $`1i,jp`$ has the entries

$$E_p(1,(i1)p+j)=\delta (i,j),$$

where $`\delta `$ is the Kronecker delta. For example, $`E_2`$ is $`[\mathrm{1\hspace{0.33em}0\hspace{0.33em}0\hspace{0.33em}1}]`$ and $`E_3`$ is $`[\mathrm{1\hspace{0.33em}0\hspace{0.33em}0\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}0\hspace{0.33em}0\hspace{0.33em}1}]`$. Let $`E_p^{}`$ be the transpose of $`E_p`$, and let us assign to $`h_k^n`$, $`1kn1`$, the $`p^n\times p^n`$ matrix $`I_{p^{nk1}}E_p^{}E_pI_{p^{k1}}`$, where $`I_m`$ is the $`m\times m`$ identity matrix with the entries $`I_m(i,j)=\delta (i,j)`$, and $``$ is the Kronecker product of matrices (see , Chapter VII.5, pp. 211-213). We assign to $`\text{1}^n`$ the matrix $`I_{p^n}`$.

Every $`n\times m`$ matrix $`A`$ whose entries are only 0 and 1 may be identified with a binary relation $`R_An\times m`$ such that $`A(i,j)=1`$ if and only if $`(i,j)R_A`$. Every binary relation may of course be drawn as a bipartite graph. Here are a few examples of such graphs for matrices we have introduced up to now, with $`p=2`$:

If $`\text{1}^n,h_1^n,\mathrm{},h_{n1}^n`$ denote the 0-1 matrices we have assigned to these expressions, then it can be verified that these matrices satisfy the equations (1), (2) and (3$`p`$), with multiplication being matrix multiplication (see , §§17-20). If $`\text{1}^n,h_1^n,\mathrm{},h_{n1}^n`$ denote the corresponding binary relations, then for multiplication being composition of binary relations the equations (1), (2) and (3) are satisfied. In both cases we also have the monoid equations.

Composition of binary relations is easy to read from bipartite graphs. Here is an example:

By so composing binary relations we can assign to every element of $`𝒥_n`$ a binary relation, and then from this binary relation we can recover the 0-1 matrix assigned to our element of $`𝒥_n`$. For this to provide an isomorphic representation of $`𝒯_n`$ in matrices it is sufficient (and also necessary) to establish that the list of 0-1 matrices assigned to the elements of $`𝒥_n`$ is linearly independent. The remainder of this paper is devoted to establishing this fact.

## 4 Diagonals

We may use $`\text{1}^n,h_1^n,\mathrm{},h_{n1}^n`$, and expressions obtained from these by multiplying, to denote either the 0-1 matrices assigned to these expressions in the previous section, or the corresponding binary relations, or the corresponding bipartite graphs. So we may speak of the 1 entries of $`h_i^n`$, conceived as a matrix, which correspond to the ordered pairs of $`h_i^n`$, conceived as a binary relation, which correspond to the edges of $`h_i^n`$, conceived as a bipartite graph. When formulating our results we will stick, however, to the terminology of binary relations (but it helps intuition if examples of such relations are drawn as bipartite graphs).

For $`n,k0`$, let $`A_k^n`$ be the set $`\{k+1,\mathrm{},k+p^n\}\text{N}`$, and for $`1qp`$ let

$$S_qA_k^{n+1}=\{(q1)p^n+jjA_k^n\}.$$

For example, if $`p=2`$, then $`A_0^2=\{1,2,3,4\}`$, $`A_0^1=\{1,2\}`$, $`A_2^1=\{3,4\}`$, $`A_2^0=\{3\}`$, $`S_1A_0^2=A_0^1`$, $`S_2A_0^2=A_2^1`$ and $`S_1S_2A_0^2=S_1(S_2A_0^2)=A_2^0`$. Words in the alphabet $`\{S_q1qp\}`$ will be called $`S`$-words.

For $`W`$ a word, let $`W^0`$ be the empty word, and let $`W^{n+1}`$ be $`W^nW`$. Let $`|W|`$ be the length of the word $`W`$. We will use the letters $`W,V,U,W_1,\mathrm{},W^{},\mathrm{}`$ for $`S`$-words.

For $`1q_1,q_2p`$, and $`q_1q_2`$, let $`L`$ be short for $`S_{q_1}`$ and $`R`$ for $`S_{q_2}`$. (When $`p=2`$, in bipartite graphs $`L=S_1`$ is interpreted as “left”, and $`R=S_2`$ as “right”.) Then the following is an easy consequence of definitions.
Remark 1.If $`n2`$ and $`|V|=n2`$, then the pairs $`(VLLA_0^n,VRRA_0^n)`$ and $`(VRRA_0^n,VLLA_0^n)`$ are in the binary relation $`h_{n1}^n`$.

From now on, we will abbreviate $`(W_1A_0^n,W_2A_0^n)`$ to $`(W_1,W_2)`$, omitting $`A_0^n`$, and, moreover, we will take for granted that such pairs belong to a binary relation, without specifying it explicitly. Then the following remark generalizes Remark 1.
Remark 2.If $`|V|0`$, $`|W|0`$, $`i=1+|V|`$ and $`n=i+1+|W|`$, then $`(VLLW,VRRW)`$ and $`(VRRW,VLLW)`$ are in $`h_i^n`$.

Next we establish the following lemma (due to the second author), which is fundamental for our proof.
Lemma 3.If $`|V|0`$, $`|W|0`$, $`k0`$, $`i=k+|V|1`$ and $`n=i+1+|W|`$, then in $`h_{i,k}^n`$ we have pairs of the following forms:
(even) for $`k=2l`$,

$$\begin{array}{cc}\text{(I }\text{even}\text{)}\hfill & (V(RL)^lLW,VL(LR)^lW)\text{,}\hfill \\ \text{(II }\text{even}\text{)}\hfill & (V(LR)^lRW,VR(RL)^lW)\text{;}\hfill \end{array}$$

(odd) for $`k=2l+1`$,

$$\begin{array}{cc}\text{(I }\text{odd}\text{)}\hfill & (V(LR)^lLLW,VRR(LR)^lW)\text{,}\hfill \\ \text{(II }\text{odd}\text{)}\hfill & (V(RL)^lRRW,VLL(RL)^lW)\text{.}\hfill \end{array}$$

(Note that (II even) is obtained from (I even) by interchanging $`L`$ and $`R`$, and analogously for (II odd) and (I odd). Note also that in (even) the $`S`$-words in between $`V`$ and $`W`$ on the right-hand sides are obtained from the corresponding words on the left-hand sides by reading them in reverse order, while in (odd) we have to read them in reverse order and interchange $`L`$ and $`R`$.)
Proof of Lemma 3. If $`k=0`$, then in $`h_{i,0}^n`$, which is $`\text{1}^n`$, we have $`(VLW,VLW)`$ and $`(VRW,VRW)`$. For $`1ki`$ we proceed by induction on $`k`$ in $`h_{i,k}^n`$. In the basis, when $`k=1`$, we apply Remark 2. For the induction step, when $`2ki`$, we have

$$h_{i,k}^n=h_{i,k1}^nh_{ik+1}^n,$$

and suppose the Lemma holds for $`h_{i,k1}^n`$.

If $`k1=2l+1`$, then for $`|V|=ik`$ and $`|W|=ni1`$ from (I odd) of the induction hypothesis we obtain that

$$(VR(LR)^lLLW,VRRR(LR)^lW)$$

is in $`h_{i,k1}^n`$, and from Remark 2 we obtain that

$$(VRRR(LR)^lW,VLLR(LR)^lW)$$

is in $`h_{ik+1}^n`$, which yields (I even) for $`h_{i,k}^n`$. We obtain analogously (II even) for $`h_{i,k}^n`$ from (II odd) for $`h_{i,k1}^n`$ and Remark 2 (we just interchange $`L`$ and $`R`$).

If $`k1=2l+2`$, then for $`|V|=ik`$ and $`|W|=ni1`$ from (I even) of the induction hypothesis we obtain that

$$(VL(RL)^{l+1}LW,VLL(LR)^{l+1}W)$$

is in $`h_{i,k1}^n`$, and from Remark 2 we obtain that

$$(VLL(LR)^{l+1}W,RR(LR)^{l+1}W)$$

is in $`h_{ik+1}^n`$, which yields (I odd) for $`h_{i,k}^n`$. We obtain analogously (II odd) for $`h_{i,k}^n`$ from (II even) for $`h_{i,k1}^n`$ and Remark 2. $`\mathrm{}`$

The pairs mentioned in Lemma 3 will be called the diagonals of the block $`h_{i,k}^n`$ for $`1ki`$. As particular cases of these diagonals we have the pairs mentioned in Remarks 1 and 2. For $`k=0`$, where $`h_{i,0}^n`$ is $`\text{1}^n`$, all the pairs of $`\text{1}^n`$ are the diagonals of $`\text{1}^n`$.

Let us say that a nonintersecting n-diagram $`D`$ of the kind introduced by Kauffman, and considered in  and , can be put into the pair $`(W_1,W_2)`$ when the threads of $`D`$ join the same $`S`$-symbols in the $`S`$-word $`W_1W_2`$. Then if $`i=k>0`$ and $`n=i+1`$, one can put into the diagonals of $`h_{i,k}^n`$ just the n-diagram corresponding to this block. An analogous property holds for the “diagonals” of other matrices, which will be defined below.

We can establish the following.
Remark 4.If $`|V|0`$, $`|W|0`$, $`k2`$, $`i=k+|V|`$, $`n=i+1+|W|`$, $`1k^{}<k`$ and $`|U|=n`$, then in $`h_{i,k^{}}^n`$ we don’t have any pair of the following forms:
(odd) for $`k=2l+1`$, $`(U,VRR(LR)^lW)`$ and $`(U,VLL(RL)^lW)`$,
(even) for $`k=2l+2`$, $`(U,VL(LR)^{l+1}W)`$ and $`(U,VR(RL)^{l+1}W)`$.
This is an easy consequence of the definition of $`h_{i,k}^n`$. It suffices to look at the right members of the diagonals of $`h_{i,k}^n`$, which cannot occur as right members of any pair of $`h_{i,k^{}}^n`$ for $`1k^{}<k`$.
Remark 5.If $`|V|0`$, $`|W|0`$, $`k0`$, $`i=k+|V|1`$, $`n=i+1+|W|`$, $`|U|=n`$ and $`|U^{}|=i`$, then
(even) for $`k=2l`$,
(I even) every pair $`(U,VL(LR)^lW)`$ of $`h_{i,k}^n`$ is either a diagonal of $`h_{i,k}^n`$ or

$`U`$ is of the form $`U^{}S_qW`$ for a symbol $`S_q`$ different from $`L`$;
(II even) every pair $`(U,VR(RL)^lW)`$ of $`h_{i,k}^n`$ is either a diagonal of $`h_{i,k}^n`$ or

$`U`$ is of the form $`U^{}S_qW`$ for a symbol $`S_q`$ different from $`R`$;
(odd) for $`k=2l+1`$,
(I odd) every pair $`(U,VRR(LR)^lW)`$ of $`h_{i,k}^n`$ is either a diagonal of $`h_{i,k}^n`$ or

$`U`$ is of the form $`U^{}S_qW`$ for a symbol $`S_q`$ different from $`L`$;
(II odd) every pair $`(U,VLL(RL)^lW)`$ of $`h_{i,k}^n`$ is either a diagonal of $`h_{i,k}^n`$ or

$`U`$ is of the form $`U^{}S_qW`$ for a symbol $`S_q`$ different from $`R`$.
This too is an easy consequence of the definition of $`h_{i,k}^n`$.

Let $`t`$ be in the normal form

$$h_{1,k_1}^n\mathrm{}h_{n1,k_{n1}}^n$$

of Section 2, and let $`\mathrm{max}(t)=\mathrm{max}\{ik_i>0\}`$. If for every $`i`$ we have $`k_i=0`$, then $`\mathrm{max}(t)`$ is undefined, and $`t`$ is equal to $`\text{1}^n`$.
Remark 6.If $`\mathrm{max}(t)=i`$ and $`n=i+1+w`$, then for every $`W`$ with $`|W|=w`$ there are $`V`$ and $`U`$ with $`|V|=|U|=i+1`$ such that in $`t`$ we have $`(VW,UW)`$, and there are no other pairs in $`t`$ except those of such a form.
This is an easy consequence of the definition of $`h_i^n`$.

If $`t`$ is in normal form as above, then we define the diagonals of $`t`$ inductively in terms of the diagonals of blocks. Formally, we define the diagonals of $`h_{1,k_1}^n\mathrm{}h_{i,k_i}^n`$ for every $`i`$ such that $`1in1`$. The diagonals of $`h_{1,k_1}^n`$ are the diagonals of this block. If $`(W_1,W_2)`$ is a diagonal of $`h_{1,k_1}^n\mathrm{}h_{j,k_j}^n`$ with $`j<n1`$ and $`(W_2,W_3)`$ is a diagonal of $`h_{j+1,k_{j+1}}^n`$, then $`(W_1,W_3)`$ is a diagonal of $`h_{1,k_1}^n\mathrm{}h_{j+1,k_{j+1}}^n`$.
Remark 7.If $`\mathrm{max}(t)=i`$ and $`n=i+1+|W|`$ for $`|W|0`$, then every diagonal of $`t`$ is of the form $`(VLW,URW)`$ or of the form $`(VRW,ULW)`$ for some $`V`$ and $`U`$ such that $`|V|=|U|=i`$.

From the definition above it is not clear that every $`t`$ has diagonals. If $`t`$ is equal to $`\text{1}^n`$, then it certainly has diagonals. For the rest we have the following lemma.
Lemma 8.Consider $`h_{1,k_1}^n\mathrm{}h_{n1,k_{n1}}^n`$ in normal form. If $`t`$ is $`h_{1,k_1}^n\mathrm{}h_{j,k_j}^n`$ for $`jn1`$, $`|V|0`$, $`|W|0`$, $`\mathrm{max}(t)=i=|V|+k_i1`$, $`n=i+1+|W|`$ and $`|U|=n`$, then, for some $`V`$ and every $`W`$, in $`t`$ we have diagonals of the following forms:
(even) for $`k_i=2l`$,

$$\begin{array}{cc}\text{(I }\text{even}\text{)}\hfill & (U,VL(LR)^lW)\text{,}\hfill \\ \text{(II }\text{even}\text{)}\hfill & (U,VR(RL)^lW)\text{;}\hfill \end{array}$$

(odd) for $`k_i=2l+1`$,

$$\begin{array}{cc}\text{(I }\text{odd}\text{)}\hfill & (U,VRR(LR)^lW)\text{,}\hfill \\ \text{(II }\text{odd}\text{)}\hfill & (U,VLL(RL)^lW)\text{.}\hfill \end{array}$$

Proof. We proceed by induction on $`i`$. For $`i=1`$ we apply Lemma 3. For $`i2`$ let $`t`$ be $`t^{}h_{i,k_i}^n\mathrm{}h_{j,k_j}^n`$, and suppose the lemma holds for $`t^{}`$. If $`t^{}`$ is equal to $`\text{1}^n`$, it suffices to apply Lemma 3 to $`h_{i,k_i}^n`$. Otherwise, let $`\mathrm{max}(t^{})=i^{}`$. We have to consider four cases.
(even-even) If $`k_i^{}=2l_i^{}`$ and $`k_i=2l_i`$, then, by the induction hypothesis, for $`|V^{}|=i^{}k_i^{}`$ and $`|W^{}|=ni^{}1`$ in $`t^{}`$ we have

$$\begin{array}{cc}(\text{I}^{}\text{even})\hfill & (U,V^{}L(LR)^{l_i^{}}W^{}),\hfill \\ (\text{II}^{}\text{even})\hfill & (U,V^{}R(RL)^{l_i^{}}W^{}),\hfill \end{array}$$

and by Lemma 3 for $`|V|=ik_i`$ and $`|W|=ni1`$ in $`h_{i,k_i}^n`$ we have

$$\begin{array}{cc}(\text{I}\text{even})\hfill & (V(RL)^{l_i}LW,VL(LR)^{l_i}W),\hfill \\ (\text{II}\text{even})\hfill & (V(LR)^{l_i}RW,VR(RL)^{l_i}W).\hfill \end{array}$$

Since $`t`$ is in normal form, we have $`i^{}<i`$ and $`i^{}k_i^{}<ik_i`$, and so we have $`|W|<|W^{}|`$ and $`|V^{}|<|V|`$. Then we can replace $`W^{}`$ in $`(\text{I}^{}\text{even})`$ and $`(\text{II}^{}\text{even})`$ by $`W^{\prime \prime }LW`$ and $`W^{\prime \prime }RW`$ for $`|W^{\prime \prime }|0`$, and we can replace $`V`$ in (I even) and (II even) by $`V^{}LV^{\prime \prime }`$ and $`V^{}RV^{\prime \prime }`$ for $`|V^{\prime \prime }|0`$. We have $`|V^{\prime \prime }|=|W^{\prime \prime }|+2(l_i^{}l_i)`$, and so $`|V^{\prime \prime }|`$ is even if and only if $`|W^{\prime \prime }|`$ is even. If both $`|V^{\prime \prime }|`$ and $`|W^{\prime \prime }|`$ are even, then we make $`(LR)^{l_i^{}}W^{\prime \prime }`$ equal to $`V^{\prime \prime }(LR)^{l_i}`$, and $`(RL)^{l_i^{}}W^{\prime \prime }`$ equal to $`V^{\prime \prime }(RL)^{l_i}`$. If both $`|V^{\prime \prime }|`$ and $`|W^{\prime \prime }|`$ are odd, then we make $`(LR)^{l_i^{}}W^{\prime \prime }`$ equal to $`V^{\prime \prime }(RL)^{l_i}`$, and $`(RL)^{l_i^{}}W^{\prime \prime }`$ equal to $`V^{\prime \prime }(LR)^{l_i}`$.

We proceed analogously in the remaining three cases, where $`k_i^{}`$ is even and $`k_i`$ odd, where $`k_i^{}`$ is odd and $`k_i`$ even, and, finally, where they are both odd. $`\mathrm{}`$

## 5 Linear independence for $`𝒥_n`$

We are now ready to prove the following result, which guarantees linear independence for the representation of $`𝒥_n`$.
Theorem.If $`t^{}`$ precedes $`t`$ in the lexicographical order, then the diagonals of $`t`$ are not pairs of $`t^{}`$.
Proof. Let $`t`$ be $`h_{1,k_1}^n\mathrm{}h_{n1,k_{n1}}^n`$, let $`t^{}`$ be $`h_{1,k_1^{}}^n\mathrm{}h_{n1,k_{n1}^{}}^n`$, and let $`t^{}`$ precede $`t`$ in the lexicographical order. Then there is an $`i`$ such that $`k_i^{}<k_i`$ and for every $`j>i`$ we have $`k_j^{}=k_j`$. We make an induction on the number $`\nu `$ of $`j`$’s such that $`j>i`$ and $`k_j^{}=k_j>0`$.

If $`\nu `$ is 0, then $`\mathrm{max}(t)=i`$. If $`t^{}`$ is $`\text{1}^n`$, then the Theorem follows by Remark 7. If $`t^{}`$ is not $`\text{1}^n`$, then let $`\mathrm{max}(t^{})=i^{}`$. Hence, $`k_i`$ and $`k_i^{}`$ are both greater than 0, and we must have $`i^{}i`$. Then two cases are possible.
(1) Suppose $`i^{}<i`$; that is, $`0=k_i^{}<k_i`$. Then, by Remark 6, every pair of $`t^{}`$ must be of the form $`(V^{}W^{},U^{}W^{})`$ for $`|V^{}|=|U^{}|=i^{}+1`$ and $`|W^{}|=ni^{}1`$. On the other hand, by Remark 7 every diagonal of $`t`$ is of the form $`(VLW,URW)`$ or $`(VRW,ULW)`$ for $`|V|=|U|=i`$ and $`|W|=ni1`$. Since $`|W|<|W^{}|`$, every diagonal of $`t`$ is not a pair of $`t^{}`$.
(2) Suppose $`i^{}=i`$; that is, $`0<k_i^{}<k_i`$. By Lemma 3, the second member of every diagonal of $`t`$ is of one of the following forms: $`VL(LR)^lW`$ and $`VR(RL)^lW`$ for $`k_i=2l`$, $`VRR(LR)^lW`$ and $`VLL(RL)^lW`$ for $`k_i=2l+1`$, where $`|V|=ik_i`$ and $`|W|=ni1`$. On the other hand, by Remark 4, no pair of $`t^{}`$ can have such second members.

Suppose now for the induction step that $`\nu >0`$. If $`\mathrm{max}(t)=j`$ and $`\mathrm{max}(t^{})=j^{}`$, then $`j=j^{}2`$ and $`k_j=k_j^{}`$. Every diagonal $`(W_1,W_3)`$ of $`t`$ can be obtained by composing a diagonal $`(W_1,W_2)`$ of $`h_{1,k_1}^n\mathrm{}h_{j1,k_{j1}}^n`$ and a diagonal $`(W_2,W_3)`$ of $`h_{j,k_j}^n`$. By the induction hypothesis, in $`h_{1,k_1^{}}^n\mathrm{}h_{j1,k_{j1}^{}}^nh_{j,0}^n\mathrm{}h_{n1,0}^n`$ we don’t have the pair $`(W_1,W_2)`$, which is a diagonal of $`h_{1,k_1}^n\mathrm{}h_{j1,k_{j1}}^nh_{j,0}^n\mathrm{}h_{n1,0}^n`$, since the former expression precedes the latter in the lexicographical order.

Suppose $`(W_1,W_3)`$ is in $`t^{}`$. Then for some $`W_4`$ we need to have $`(W_1,W_4)`$ in $`h_{1,k_1^{}}^n\mathrm{}h_{j1,k_{j1}^{}}^n`$ and $`(W_4,W_3)`$ in $`h_{j,k_j}^n`$. It is excluded that $`W_4`$ is $`W_2`$, as we remarked above.

Suppose $`k_j=2l+1`$, and suppose $`(W_2,W_3)`$ is $`(V(LR)^lLLW,VRR(LR)^lW)`$ for $`|W|=nj1`$. By Lemma 3, this is one possible case. Then, by Remark 5, we must have that $`W_4`$ is $`W_4^{}S_qW`$ for $`S_q`$ different from $`L`$, since $`W_4`$ cannot be $`W_2`$. By Remark 7, we have that $`W_1`$ must be of the form $`W_1^{}LW`$. By Remark 6, we cannot have $`(W_1^{}LW,W_4^{}S_qW)`$ in $`h_{1,k_1^{}}^n\mathrm{}h_{j1,k_{j1}^{}}^n`$. So $`(W_1,W_3)`$ is not in $`t^{}`$. We proceed analogously in other possible cases. $`\mathrm{}`$

As a corollary of Lemma 8 and of the Theorem, we obtain that every matrix in $`𝒥_n`$ has an entry with 1 (one of its “diagonals”) where all the matrices preceding it in the lexicographical order have 0. This implies that the matrices assigned to the elements of $`𝒥_n`$ make a linearly independent list, and hence Brauer’s representation of Temperley-Lieb algebras in matrices, which we have considered here, is faithful.


text:
# XX QCD Coherence Effects in Low-x DIS and Drell-Yan with Nuclear Targets

## I ABSTRACT

We investigate QCD coherence effects in deep inelastic scattering (DIS) off nuclei and in Drell-Yan (DY) dilepton production in proton-nucleus collisions within the light-cone color-dipole approach. The physical mechanisms underlying the nuclear effects become very transparent in this approach and are explained in some detail. We present numerical calculations of nuclear shadowing in DIS and DY and compare to data. Nuclear effects in the DY transverse momentum distribution are calculated as well. The dipole approach is the only known way to calculate the Cronin effect without introducing additional parameters for nuclear targets.

Keywords: Deep Inelastic Scattering, Drell-Yan Process, QCD, Nuclear Effects, Coherence Effects, Multiple Scattering.

## II 1. INTRODUCTION

The use of nuclei instead of protons in high energy scattering experiments, like deep inelastic scattering, provides unique possibilities to study the space-time development of strongly interacting systems. In experiments with proton targets the products of the scattering process can only be observed in a detector which is separated from the reaction point by a macroscopic distance. In contrast to this, the nuclear medium can serve as a detector located directly at the place where the microscopic interaction happens. As a consequence, with nuclei one can study coherence effects in QCD which are not accessible in DIS off protons nor in proton-proton scattering.

At high energies, nuclear scattering is governed by coherence effects which are most easily understood in the target rest frame. In the rest frame, DIS looks like pair creation from a virtual photon, see Fig. 1. Long before the target, the virtual photon splits into a $`q\overline{q}`$-pair. The lifetime $`l_c`$ of the fluctuation, which is also called coherence length, can be estimated with help of the uncertainty relation to be of order $`1/m_Nx_{Bj}`$, where $`x_{Bj}`$ is Bjorken-x and $`m_N1`$ GeV is the nucleon mass. The coherence length can become much greater than the nuclear radius at low $`x_{Bj}`$. Multiple scattering within the lifetime of the $`q\overline{q}`$ fluctuation leads to the pronounced coherence effects observed in experiment.

The most prominent example for a coherent interaction of more than one nucleon is the phenomenon of nuclear shadowing, i.e. the suppression of the nuclear structure function $`F_2^A`$ with respect to the proton structure function $`F_2^p`$ at low $`x_{Bj}\text{ }<0.1`$, $`F_2^A(x_{Bj},Q^2)/(AF_2^p(x_{Bj},Q^2))`$. Here $`Q^2`$ is the virtuality of the photon. Shadowing in low $`x_{Bj}`$ DIS and at high photon virtualities is experimentally well studied by NMC NMC . The analogous effect for DY at low $`x_2`$ was discovered by E772 E772 .

What is the mechanism behind this suppression? If the coherence length is very long, as indicated in Fig. 1, the $`q\overline{q}`$-dipole undergoes multiple scatterings inside the nucleus. The physics of shadowing in DIS is most easily understood in a representation, in which the pair has a definite transverse size $`\rho `$. As a result of color transparency zkl ; bbgg ; bm , small pairs interact with a small cross section, while large pairs interact with a large cross section. The term ”shadowing” can be taken literally in the target rest frame. Large pairs are absorbed by the nucleons at the surface which cast a shadow on the inner nucleons. The small pairs are not shadowed. They have equal chances to interact with any of the nucleons in the nucleus. From these simple arguments, one can already understand the two necessary conditions for shadowing. First, the hadronic fluctuation of the virtual photon has to interact with a large cross section and second, the coherence length has to be long enough to allow for multiple scattering.

While the target rest frame picture of DIS is very popular, the light-cone approach of Kopeliovich boris ; bhq ; kst1 , which describes the DY process dy  in the target rest frame, is less known. In the light-cone approach, DY dilepton production in the rest frame of the target appears as bremsstrahlung, see Fig. 2. A quark from the projectile scatters off the target and radiates a virtual photon. This photon decays into a lepton pair. Even though this mechanism looks quite different from the usual annihilation picture, it yields numerical values very similar to the next-to-leading order parton model new . Note that while cross sections are Lorentz invariant, the partonic interpretation of a hard scattering process is not. That is why a process that looks like quark-antiquark annihilation in the infinite momentum frame can appear as bremsstrahlung in the target rest frame.

Remarkably, the DY cross section can be expressed in terms of the same dipole cross section that appears in DIS. In DY, the dipole cross section arises from the interference of the two graphs in Fig. 2. On a nuclear target, the quark will of course scatter several times. The effect of multiple scattering on bremsstrahlung is well known in QED as the Landau-Pomeranchuk-Migdal (LPM) effect LPM . The LPM effect leads to a reduction of the cross section due to destructive interferences between the amplitudes for radiation off different points along the quark’s trajectory. In the case of DY, this suppression is observed in experiment as nuclear shadowing E772 .

## III 2. SHADOWING IN DIS

First, we discuss the case of shadowing in deep inelastic scattering off nuclei. Like shadowing in hadron-nucleus collisions glaubergribov , also shadowing in DIS is intimately related to diffraction. The close connection between shadowing and diffraction becomes most transparent in the formula derived by Karmanov and Kondratyuk KKK . In the double scattering approximation, the shadowing correction can be related to the diffraction dissociation spectrum, integrated over the mass,

$`\sigma ^{\gamma ^{}A}`$ $``$ $`A\sigma ^{\gamma ^{}p}`$
$``$ $`4\pi {\displaystyle 𝑑M_X^2\frac{d\sigma (\gamma ^{}NXN)}{dM_X^2dt}}|_{t0}`$ (1)
$`\times `$ $`{\displaystyle d^2bF_A^2(l_c,b)}.`$
Here

$$F_A^2(l_c,b)=\left|\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑z\rho _A(b,z)\mathrm{e}^{\mathrm{i}z/l_c}\right|^2$$
(2)
is the formfactor of the nucleus, which depends on the coherence length

$$l_c=\frac{2\nu }{Q^2+M_X^2},$$
(3)
$`\nu `$ is the energy of the $`\gamma ^{}`$ in the target rest frame and $`M_X`$ is the mass of the diffractively excited state. The coherence length can be estimated from the uncertainty relation and is the lifetime of the diffractively excited state. If $`l_c0`$, the shadowing correction in Eq. (1) vanishes and one is left with the single scattering term $`A\sigma ^{\gamma ^{}p}`$.

Note that Eq. (1) is valid only in double scattering approximation. For heavy nuclei, however, higher order scattering terms will become important. These can be calculated, if one knows the eigenstates of the interaction. Fortunately, the eigenstates of the $`T`$ matrix (restricted to diffractive processes) were identified a long time ago in QCD zkl ; mp  as partonic configurations with fixed transverse separations in impact parameter space. For DIS, the lowest eigenstate is the $`q\overline{q}`$ Fock component of the photon. The total $`\gamma ^{}`$-proton cross section is easily calculated, if one knows the cross section $`\sigma _{q\overline{q}}(\rho )`$ for scattering a $`q\overline{q}`$-dipole of transverse size $`\rho `$ off a proton,

$$\sigma ^{\gamma ^{}p}=𝑑\alpha d^2\rho \left|\mathrm{\Psi }_{q\overline{q}}(\alpha ,\rho )\right|^2\sigma _{q\overline{q}}(\rho ).$$
(4)
The light-cone wavefunction $`\mathrm{\Psi }_{q\overline{q}}(\alpha ,\rho )`$ describes the splitting of the virtual photon into the $`q\overline{q}`$-pair and is calculable in perturbation theory, see e.g. nz . Here, $`\alpha `$ is the longitudinal momentum fraction carried by the quark in Fig. 1. The dipole cross section is governed by nonperturbative effects and cannot be calculated from first principles. We use the phenomenological parameterization from thesis  for our calculation of shadowing in DIS, which is fitted to HERA data on the proton structure function. Note that higher Fock-states of the photon, containing gluons, lead to an energy dependence of $`\sigma _{q\overline{q}}`$, which we do not write out explicitly.

The diffractive cross section can also be expressed in terms of $`\sigma _{q\overline{q}}`$. Since the cross section for diffraction is proportional to the square of the $`T`$-matrix element, $`|\gamma ^{}|T|X|^2`$, the dipole cross section also enters squared,

$$𝑑M_X^2\frac{d\sigma (\gamma ^{}NXN)}{dM_X^2dt}|_{t0}=\frac{\sigma _{q\overline{q}}^2(\rho )}{16\pi },$$
(5)
where the brackets $`\mathrm{}`$ indicate averaging over the light-cone wavefunctions like in Eq. (4). We point out that in order to reproduce the correct behavior of the diffractive cross section at large $`M_X`$, one has to include at least the $`q\overline{q}G`$ Fock-state of the $`\gamma ^{}`$. This correction is, however, of minor importance in the region where shadowing data are available.

If one attempts to calculate shadowing from Eq. (1) with help of Eq. (5), one faces the problem that the nuclear form factor, Eq. (2), depends on the mass $`M_X`$ of the diffractively produced state, which is undefined in impact parameter representation. Only in the limit $`l_cR_A`$, where $`R_A`$ is the nuclear radius, it is possible to resum the entire multiple scattering series in an eikonal-formula

$$\sigma ^{\gamma ^{}A}=2d^2b\left(1\mathrm{exp}\left(\frac{\sigma _{q\overline{q}}(\rho )}{2}T(b)\right)\right).$$
(6)
The nuclear thickness function $`T(b)=_{\mathrm{}}^{\mathrm{}}𝑑z\rho _A(b,z)`$ is the integral of nuclear density over longitudinal coordinate $`z`$ and depends on the impact parameter $`b`$. The condition $`l_cR_A`$ insures that the $`\rho `$ does not vary during propagation through the nucleus (Lorentz time dilation) and one can apply the eikonal approximation.

The condition $`l_cR_A`$ is however not fulfilled in experiment. For the case $`l_cR_A`$, one has to take the variation of $`\rho `$ during propagation of the $`q\overline{q}`$ fluctuation through the nucleus into account, see Fig. 3. A widely used recipe is to replace $`M_X^2Q^2`$, so that $`l_c1/(2m_Nx_{Bj})`$ and one can apply the double scattering approximation. This recipe was, however, disfavored during our investigation krt2 . Moreover, there is no simple recipe to include a finite $`l_c`$ into higher order scattering terms.

In krt1  a Green function technique was developed that provides the correct quantum-mechanical treatment of a finite coherence length in all multiple scattering terms. Like in Eq. (1) the total cross section is represented in the form

$$\sigma ^{\gamma ^{}A}=A\sigma ^{\gamma ^{}p}\mathrm{\Delta }\sigma ,$$
(7)
where $`\mathrm{\Delta }\sigma `$ is the shadowing correction,

$`\mathrm{\Delta }\sigma `$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Re}{\displaystyle d^2b\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑z_1\rho _A(b,z_1)\underset{z_1}{\overset{\mathrm{}}{}}𝑑z_2\rho _A(b,z_2)}`$ (8)
$`\times `$ $`{\displaystyle _0^1}𝑑\alpha {\displaystyle d^2\rho _2\mathrm{\Psi }_{q\overline{q}}^{}(\stackrel{}{\rho }_2,\alpha )\sigma _{q\overline{q}}(\rho _2)}`$
$`\times `$ $`A(\stackrel{}{\rho }_2,z_1,z_2,\alpha ),`$
with

$`A(\stackrel{}{\rho }_2,z_1,z_2,\alpha )`$ $`=`$ $`{\displaystyle d^2\rho _1W(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)}`$
$`\times `$ $`\mathrm{e}^{\mathrm{i}q_L^{min}(z_2z_1)}\sigma _{q\overline{q}}(\rho _1)\mathrm{\Psi }_{q\overline{q}}(\stackrel{}{\rho }_1,\alpha ).`$
Here,

$$q_L^{min}=\frac{1}{l_c^{max}}=\frac{Q^2\alpha (1\alpha )+m_f^2}{2\nu \alpha (1\alpha )}$$
(10)
is the minimal longitudinal momentum transfer when the photon splits into the $`q\overline{q}`$ dipole ($`m_f`$ is the quark mass).

The shadowing term in Eq. (7) is illustrated in Fig. 3. At the point $`z_1`$ the photon diffractively produces the $`q\overline{q}`$ pair ($`\gamma ^{}Nq\overline{q}N`$) with a transverse separation $`\stackrel{}{\rho }_1`$. The pair propagates through the nucleus along arbitrarily curved trajectories, which are summed over, and arrives at the point $`z_2`$ with a separation $`\stackrel{}{\rho }_2`$. The initial and the final separations are controlled by the light-cone wavefunction $`\mathrm{\Psi }_{q\overline{q}}(\stackrel{}{\rho },\alpha )`$. While passing the nucleus the $`q\overline{q}`$ pair interacts with bound nucleons via the cross section $`\sigma _{q\overline{q}}(\rho )`$ which depends on the local separation $`\stackrel{}{\rho }`$. The Green function $`W(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)`$ describes the propagation of the pair from $`z_1`$ to $`z_2`$, see Eq. (III), including all multiple rescatterings and a finite coherence length. Note the diffraction dissociation ($`DD`$) amplitude,

$$f_{DD}(\gamma ^{}q\overline{q})=\mathrm{i}\mathrm{\Psi }_{q\overline{q}}(\stackrel{}{\rho }_1,\alpha )\sigma _{q\overline{q}}(\rho _1)$$
(11)
, in Eq. (III). At the position $`z_2`$, the result of the propagation is again projected onto the diffraction dissociation amplitude, Eq. (8). The Green function includes that part of the phase shift between the initial and the final photons which is due to transverse motion of the quarks, while the longitudinal motion is included in Eq. (III) via the exponential.

The Green function $`W(\stackrel{}{\rho }_2,z_2;\stackrel{}{\rho }_1,z_1)`$ in Eq. (III) satisfies the two dimensional Schrödinger equation,

$`\mathrm{i}{\displaystyle \frac{W(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)}{z_2}}=`$ (12)
$``$ $`{\displaystyle \frac{\mathrm{\Delta }(\rho _2)}{2\nu \alpha (1\alpha )}}W(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)`$
$``$ $`{\displaystyle \frac{\mathrm{i}}{2}}\sigma (\rho _2)\rho _A(b,z_2)W(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)`$
with the boundary condition $`W(\stackrel{}{\rho }_2,z_1|\stackrel{}{\rho }_1,z_1)=\delta ^{(2)}(\stackrel{}{\rho }_2\stackrel{}{\rho }_1)`$. The Laplacian $`\mathrm{\Delta }(\rho _2)`$ acts on the coordinate $`\stackrel{}{\rho }_2`$. The kinetic term $`\mathrm{\Delta }/[2\nu \alpha (1\alpha )]`$ in this Schrödinger equation takes care of the varying effective mass of the $`q\overline{q}`$ pair and provides the proper phase shift. The role of time is played by the longitudinal coordinate $`z_2`$. The imaginary part of the optical potential describes the rescattering.

The Green function method contains the eikonal approximation Eq. (6) and the Karmanov-Kondratyuk formula Eq. (1) as limiting cases. In order to obtain the eikonal approximation, one has to take the limit $`\nu \mathrm{}`$. In this case, the kinetic energy term in Eq. (12) can be neglected and with $`q_L^{min}0`$ one arrives after a short calculation at Eq. (6).

One can also recover the Karmanov-Kondratyuk formula, when one neglects the imaginary potential in Eq. (12). Then $`W`$ becomes the Green function of a free motion.

Calculations in the Green function approach are compared to NMC data in Fig. 4. Note that the leading twist contribution to shadowing is due to large dipole sizes, where nonperturbative effects, like an interaction between the $`q`$ and the $`\overline{q}`$, might become important. Therefore, the solid curve is calculated following kst2 , including such an interaction. This interaction modifies the light-cone wavefunction kst2 . The dashed curve is calculated with the conventional, perturbative light-cone wavefunctions, but including a constituent quark mass. Both curves are in reasonable agreement with the data.

Note that for these data, the coherence length is of order of the nuclear radius or smaller. Indeed, shadowing vanishes around $`x_{Bj}0.1`$, because the coherence length becomes smaller than the mean internucleon spacing. Therefore, the eikonal approximation, Eq. (6), cannot be applied for the kinematics of NMC and a correct treatment of the coherence length becomes crucial. We emphasize that the calculation in Fig. 4 does not contain any free parameters. Following the spirit of Glauber theory, all free parameters are adjusted to DIS off protons.

## IV 3. NUCLEAR EFFECTS IN THE DRELL-YAN PROCESS

In the second part of this paper, we discuss the application of the dipole approach to calculate nuclear effects in DY. The main features of the approach have already been introduced in the preceding section and the case of DY is almost completely analogous to DIS.

The cross section for radiation of a virtual photon from a quark after scattering on a proton, can be written in factorized light-cone form boris ; bhq ; kst1 ,

$$\frac{d\sigma (qp\gamma ^{}X)}{d\mathrm{ln}\alpha }=d^2\rho |\mathrm{\Psi }_{\gamma ^{}q}(\alpha ,\rho )|^2\sigma _{q\overline{q}}(\alpha \rho ),$$
(13)
with the same $`\sigma _{q\overline{q}}`$ as in DIS. Again, the light-cone wavefunction $`\mathrm{\Psi }_{\gamma ^{}q}(\alpha ,\rho )`$ is calculable in perturbation theory krt3 ; new . For DY, however, $`\rho `$ is the photon-quark transverse separation and $`\alpha `$ is the fraction of the light-cone momentum of the initial quark taken away by the photon. We use the standard notation for the kinematical variables, $`x_1x_2=x_F`$, $`\tau =M^2/s=x_1x_2`$, where $`x_F`$ is the Feynman variable, $`s`$ is the center of mass energy squared of the colliding protons and $`M`$ is the dilepton mass.

The physical interpretation of (13) is similar to the DIS case. Long before the target, the projectile quark develops fluctuations which contain virtual photons. On interaction with the target, these fluctuations can be put on mass shell and the $`\gamma ^{}`$ is freed. In order to make this happen, the interaction must distinguish between Fock states with and without $`\gamma ^{}`$, i.e. these states have to interact differently. Since only the quark interacts in both Fock components, the difference in the scattering amplitudes arises from the relative displacement of the quark in the transverse plane when it radiates a photon. The $`\gamma ^{}q`$ fluctuation has a center of gravity in the transverse plane which coincides with the impact parameter of the parent quark. One finds that the transverse separation between the quark and the center of gravity is $`\alpha \rho `$, i.e. the argument of $`\sigma _{q\overline{q}}`$. The dipole cross section arises from the interference of the two graphs in Fig. 2 and should be interpreted as freeing cross section for the virtual photon. More discussion can be found in krt3 .

The transverse momentum distribution of DY pairs can also be expressed in terms of the dipole cross section kst1 . The differential cross section is given by the Fourier integral

$`{\displaystyle \frac{d\sigma (qp\gamma ^{}X)}{d\mathrm{ln}\alpha d^2q_T}}=`$
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle d^2\rho _1d^2\rho _2\mathrm{exp}[\mathrm{i}\stackrel{}{q}_T(\stackrel{}{\rho }_1\stackrel{}{\rho }_2)]}`$
$`\times `$ $`\mathrm{\Psi }_{\gamma ^{}q}^{}(\alpha ,\stackrel{}{\rho }_1)\mathrm{\Psi }_{\gamma ^{}q}(\alpha ,\stackrel{}{\rho }_2)`$
$`\times `$ $`{\displaystyle \frac{1}{2}}\left\{\sigma _{q\overline{q}}(\alpha \rho _1)+\sigma _{q\overline{q}}(\alpha \rho _2)\sigma _{q\overline{q}}(\alpha (\stackrel{}{\rho }_1\stackrel{}{\rho }_2))\right\}.`$
After integrating this expression over the transverse momentum $`q_T`$ of the photon, one obviously recovers Eq. (13).

Nuclear effects arise in the DY process, because the different nucleons in the nucleus compete to free the $`\gamma ^{}`$. Necessary condition for this phenomenon is of course that the coherence length for DY,

$$l_c=\frac{1}{m_Nx_2}\frac{(1\alpha )M^2}{q_T^2+(1\alpha )M^2+\alpha ^2m_f^2}$$
(15)
exceeds the mean internucleon spacing of $`2`$ fm. At infinite coherence length, one can simply eikonalize the dipole cross section, similar to the DIS case, cmp. Eq. (6). For finite coherence length, the Green function technique can be developed for DY in complete analogy to DIS. We do not display the corresponding formula here, which can be found in kst1 ; thesis ; moriond . We point out that even though DIS and DY appear very similar in the dipole approach, there are no (or only very few) diffractively produced DY pairs, because Eq. (5) is not applicable to an exclusive channel (see discussion in kst2 ).

Calculations of shadowing for DY in the Green function approach are compared to E772 data E772  in Fig. 5. Like in the case of the NMC data in the previous section, the coherence length $`l_c`$ at E772 energy becomes smaller than the nuclear radius. Shadowing vanishes as $`x_2`$ approaches $`0.1`$, because the coherence length becomes smaller than the mean internucleon separation. Again, it is therefore important to have a correct description of a finite $`l_c`$. The eikonal approximation does not reproduce the vanishing shadowing toward $`x_20.1`$, because it assumes an infinite coherence length.

The curves in Fig. 5 are somewhat different from the ones in thesis ; moriond , because we used a different parameterization of the dipole cross section wuesthoff . Note that for heavy nuclei, energy loss eloss  leads to an additional suppression of the DY cross section. It is therefore important to calculate shadowing in DY without introducing additional parameter for nuclear targets. Otherwise, shadowing and energy loss cannot be disentangled.

Nuclear effects on the $`q_T`$-differential cross section calculated at RHIC and LHC energy are shown in Fig. 6. At these high cm energies, the coherence length significantly exceeds the nuclear radius, justifying the eikonal approximation. However, in the case of such long $`l_c`$, the quark has sufficient time to radiate additional gluons. Thus, higher Fock states have to be taken into account. This is accomplished by multiplying $`\sigma _{q\overline{q}}`$ with the gluon shadowing ratio $`R_G`$, which was calculated in kst2 . See krtj  for a detailed discussion.

The differential cross section is suppressed at small transverse momentum $`q_T`$ of the dilepton, where large values of $`\rho `$ dominate. This suppression vanishes at intermediate $`q_T2`$ GeV. The Cronin enhancement that one could expect in this intermediate $`q_T`$ region moriond  is suppressed due to gluon shadowing krtj . At very large transverse momentum nuclear effects vanish.

## V 4. SUMMARY

In the target rest frame and at low $`x_{Bj}`$ ($`x_2`$), the cross sections for DIS and DY can both be expressed in terms of the same cross section for scattering a color-neutral $`q\overline{q}`$-pair off a proton. The advantage of this formulation and the main motivation to develop this approach is the easy generalization to nuclear targets and the insight into the dynamical origin of nuclear effects, which appear as coherence effects due to multiple scattering.

The main nonperturbative input to all formulae is the dipole cross section, which cannot be calculated from first principles. Instead one uses a phenomenological parameterization for this quantity, which is determined to fit DIS data from HERA. Nuclear effects are then calculated without introducing any new parameters. A parameter free calculation of nuclear shadowing for DY is indispensable, if one aims at extracting the energy loss of a fast quark propagating through nuclear from DY data eloss . Furthermore, the dipole formulation is the only known way to calculate the Cronin effect in a parameter free way. This approach can also be applied to calculate the Cronin effect in hadronic collisions ce . Furthermore, a good description of NMC data on shadowing in DIS is achieved.

Finally, we mention that the dipole approach can also be applied to a variety of other processes, e.g. heavy quark, quarkonium or vector meson production bk .

Acknowledgments: This work was supported by the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. W-7405-ENG-38. The author is grateful to Jörg Hüfner, Mikkel Johnson, Boris Kopeliovich, Jen-Chieh Peng and Alexander Tarasov for valuable discussion.


text:
# Fast crack propagation by surface diffusion

## Abstract

We present a continuum theory which describes the fast growth of a crack by surface diffusion. This mechanism overcomes the usual cusp singularity by a self-consistent selection of the crack tip radius. It predicts the saturation of the steady state crack velocity appreciably below the Rayleigh speed and tip blunting. Furthermore, it includes the possibility of a tip splitting instability for high applied tensions.

One of the most challenging puzzles in nonequilibrium physics and materials science is the phenomenon of fracture. It is important for the vast field of material failure and probably also for friction processes Gerde01 . Despite its relevance even the motion of a single crack is poorly understood Hellemans98 . Experimentally, the maximum attained crack velocities are far lower than the theoretically expected Rayleigh speed Freund90 . Beyond a critical velocity, a so far unpredictable tip splitting of the crack can happen and produce strange oscillations of the crack speed Sharon99 .

The classical theories Freund90  are based on linear theory of elasticity and on an integral energy balance in the vicinity of an infinitely sharp crack tip. However, a more detailed approach based on equations of motion for the crack shape is needed to describe the intriguing spectrum of phenomena near the crack tip. Hence, the curvature of the crack tip is required as a new relevant dynamical variable which also allows to avoid stress singularities. We emphasize that, in contrast to models which describe crack propagation by bond breaking at the infinitely sharp tip, growth with a finite tip radius always requires a transport mechanism in order to preserve the shape (see Fig. 1). Recently it was proposed that this lengthscale is dynamically selected by the threshold of plastic deformations in the tip region Langer00 . Unfortunately, approaches of this type (see also Aranson00 ) require the introduction of theories of plasticity which are usually much more speculative and less verified than the ordinary linear theory of elasticity.

Here we demonstrate that linear theory of elasticity is sufficient to describe consistently crack propagation, driven by surface diffusion along the crack surfaces. Of course, in many situations plasticity is very important, but the beauty of our approach is that it predicts, in a simple and well controlled continuum theory, steady state crack growth, the tip splitting instability and also slow deformations of already existing cracks.

The idea of crack propagation by surface diffusion has previously been studied by Stevens and Dutton Stevens71 , who assumed an ad hoc crack shape which is not found by solving the full free boundary problem and requires mass transport over large scales. Therefore, their model cannot describe the usual fast crack growth.

Our basic mechanism is related to the Asaro-Tiller-Grinfeld (ATG) instability Asaro72 , which predicts a morphological instability of a uniaxially stressed solid interface due to surface diffusion. Relatively long-wave perturbations of the interface lead to a reduction of the elastic energy of the system, whereas short-wave corrugations are hampered by surface energy. In the long time behavior, deep grooves can form, producing shapes similar to cracks Yang93 ; Spencer94 . According to previous theories, which used only the static theory of elasticity, the notches propagate with increasing velocity and decreasing tip radius and collapse to a finite-time cusp singularity. Similar to crack dynamics the lack of tip radius selection becomes obvious, and already shows the close relationship between the ATG instability and crack propagation.

Usually, it is believed that surface diffusion is slow, but, surprisingly enough, it should not be ignored even in fast fracture processes. Our main idea is that surface diffusion is driven by the strong gradient of the chemical potential in the tip region. This can be a very efficient mechanism for crack propagation if the transport length is sufficiently small. Additionally, energy release and strong dissipation bring the local temperature close to the melting temperature Fuller75 ; Fineberg99 . This drastically increases the surface diffusion coefficient and makes fast crack propagation essentially independent of the outside temperature.

Steady state crack growth — We use a two-dimensional plane-strain situation with mode I loading to describe crack propagation Freund90 .

On the surfaces of the crack the normal stress $`\sigma _{nn}`$ and the shear stress $`\sigma _{n\tau }`$ vanish, whereas the tangential stress $`\sigma _{\tau \tau }`$ usually does not. The chemical potential at the interface is given by

$$\mu =\mathrm{\Omega }\left(\frac{1\nu ^2}{2E}\sigma _{\tau \tau }^2\alpha \kappa \right).$$
(1)
Here $`\alpha `$ is the surface energy, $`\kappa `$ the curvature of the interface and $`\mathrm{\Omega }`$ the atomic volume. $`E`$ and $`\nu `$ are Young’s modulus and Poisson ratio, respectively. Non-hydrostatic stresses drive a surface flux proportional to the gradients of the chemical potential along the surface; in turn the normal velocity equals the divergence of this flux due to conservation of material,

$$v_n=\frac{D}{\alpha \mathrm{\Omega }}\frac{^2\mu }{s^2}$$
(2)
where $`/s`$ denotes the tangential derivative and $`D`$ (dimension $`\mathrm{m}^4\mathrm{s}^1`$) is proportional to the surface diffusion coefficient. (It is related to the usual surface diffusion coefficient $`D_s`$ by $`D=D_s\mathrm{\Omega }^2\delta \alpha /kT`$. Here $`\delta `$ is the number of atoms per unit area of surface, $`k`$ the Boltzmann constant and $`T`$ the temperature.)

First, we are interested in steady state solutions of the equation of motion, with a crack moving in positive $`x`$-direction with velocity $`v`$ (see Fig. 1); In co-moving polar coordinates, $`x=r(\theta )\mathrm{cos}\theta ,y=r(\theta )\mathrm{sin}\theta `$, the steady state equation for the shape $`r(\theta )`$ reads after one integration of Eq. (2)

$$vr\mathrm{sin}\theta =\frac{D}{\alpha \mathrm{\Omega }}\frac{1}{\sqrt{r^2+r^2}}\frac{d\mu }{d\theta }.$$
(3)
Generally speaking, this, together with Eq. (1), is a complicated, non-linear third order equation with non-local contributions arising from the elastic fields, since $`\sigma _{\tau \tau }`$ depends on the entire shape.

In the tail region stresses decay and the shape equation is a third order linear differential equation $`Dy^{\prime \prime \prime }=vy`$ with two growing (and oscillating) and one decaying solution. Only the latter, $`y(x\mathrm{})=A\mathrm{exp}[(v/D)^{1/3}x]`$, asymptotically describes physical shapes and is allowed. Let us focus on symmetrical solutions, $`r(\theta )=r(\theta )`$, and start integration at the crack tip $`\theta =0`$. Since the physical properties, curvature and stresses, do not depend on the choice of coordinate system but only on the crack shape, we can arbitrarily chose $`r(\theta =0)=r_0`$, with the a priori unknown tip radius $`r_0=1/\kappa (0)`$. Then from symmetry and the definition of the tip curvature, $`\kappa =(r^2+2r^2rr^{\prime \prime })/(r^2+r^2)^{3/2}`$, the natural conditions $`r^{}(0)=r^{\prime \prime }(0)=0`$ arise. Integration over the upper interface $`\theta >0`$ requires the suppression of two growing exponentials at the tail, which imposes two boundary conditions. For a given external loading, these conditions can be fulfilled by a proper selection of the tip radius $`r_0`$ and growth velocity $`v`$. By this argument the situation seems to be fully described. However, as we mentioned already earlier, the use of a static theory of elasticity does not allow a selection of the tip radius. The reason is, that both contributions to the chemical potential, surface energy $`\mu _s\kappa `$ and elastic energy $`\mu _{el}\sigma ^2`$, behave as $`r_0^1`$ close to the tip: In the tip approximation, stresses behave as Freund90

$$\sigma _{ij}=\frac{K}{r^{1/2}}f_{ij}(\theta )$$
(4)
with the static stress intensity factor $`K\sigma _{\mathrm{}}L^{1/2}`$, where $`\sigma _{\mathrm{}}`$ is the applied remote stress and $`L`$ is the macroscopic length of the crack, $`Lr_0`$, which is not considered here. Instead we assume $`K`$ to be kept fixed. The universal stress distribution $`f_{ij}`$ depends only on the orientation relative to the crack Freund90 . The asymptotic distribution (4) is valid far away from the tip, $`r_0rL`$, but it gives the correct scaling of stresses also on the crack surface $`rr_0`$. Therefore, a dimensionless rescaling of all lengthscales, e.g. $`\stackrel{~}{r}=r/r_0`$, and of the growth velocity $`\stackrel{~}{v}=vr_0^3/D`$ leaves the equation of motion invariant and, thus, cannot determine the lengthscale $`r_0`$. Consequently, a steady state solution does not exist. This is the reason for the already mentioned cusp singularity of the ATG instability.

The main idea of this paper is based on the fact that a full elastodynamic description restores the selection of this lengthscale. It is known that at least for higher crack speeds the angular distribution $`f_{ij}`$ become strongly dependent on the ratio $`v/v_R`$ ($`v_R`$ is the Rayleigh speed Freund90 ). The dynamical stress intensity factor $`K_{dyn}`$ is related to the static one used here by an extra velocity dependent function $`g(v/v_R)`$, $`K_{dyn}=Kg(v/v_R)`$. The crucial observation is, that velocity appears now in two different combinations in the equation of motion, $`vr_0^3/D`$ and $`v/v_R`$. Thus, by introduction of the new parameter $`v/v_R`$, a selection of both $`v`$ and $`r_0`$ happens.

From these general arguments we conclude that fast steady state crack propagation by surface diffusion is indeed possible. However, the exact solution of the problem is technically very difficult, because it requires the solution of an elastodynamic problem for an a priori unknown crack shape. The bulk equations of elasticity

$$\frac{\sigma _{ij}}{x_j}=\rho \ddot{u}_i$$
(5)
are subject to the boundary conditions on the crack surface (surface of discontinuity) Freund90

$$\sigma _{in}+\rho \dot{u}_iv_n=0$$
(6)
with the mass density $`\rho `$ and the elastic displacement field $`u_i`$. This is just the momentum balance equation on the free surface which moves with normal velocity $`v_n`$. Finally, the expression for the chemical potential on the crack surface should also be corrected compared to the static case:

$$\mu =\mathrm{\Omega }\left(\frac{1}{2}\sigma _{ik}u_{ik}\frac{1}{2}\rho \dot{u}_i^2\alpha \kappa \right)$$
(7)
Eqs. (5)-(7) together with surface diffusion Eq. (2) and given loading configuration describe the crack propagation in our model. We note that inertial effects appear not only in the bulk equations of elasticity but also in the boundary conditions and in the expression for the chemical potential. All these effects lead to the appearance of the parameter $`(v/v_R)^2`$ in the problem compared to the quasistatic description.

The preceding equations (5)-(7) can be derived from the Lagrangian

$$=_{V(t)}\left(\frac{1}{2}\rho \dot{u}_i^2\frac{1}{2}\sigma _{ij}u_{ij}\right)𝑑V_{S(t)}\alpha 𝑑S.$$
(8)
with $`V(t)`$ being the time dependent volume of the solid and $`S(t)`$ its surface. The elastic equations and boundary conditions follow by the condition that $`𝒮=𝑑t`$ is stationary with respect to variations of the displacement $`u_i`$, and the chemical potential is related to the variation of $``$ with respect to the interface position Uwaha86 ; Unpub . Evaluation of the exact equations of motion requires extended numerics, especially in the time-dependent case.

The local crack tip model — We simplify the problem in order to make further analytical progress and to expose the general idea of our approach. It will turn out that one cannot describe all effects by this approximation and further refinement is necessary, but the main results are qualitatively very robust against changes of the model.

We mimic the tangential stress by a local description in the spirit of Eq. (4), as depending on the propagation velocity and only on the local properties of the interface. It takes both the velocity dependence of the angular distribution and the decrease of the dynamical stress intensity factor into account:

$$\sigma _{\tau \tau }=K\left[\sqrt{1(v/v_R)^2}\mathrm{cos}(\theta /2)+(v/v_R)^2\mathrm{sin}^4\theta \right]/r^{1/2}.$$
(9)
This form reflects the first order transition of the principal stress direction $`\theta =0`$ for low velocities towards $`\theta 0`$ as function of $`v/v_R`$ Freund90 . The use of more sophisticated expressions (e.g. the singular dynamical field in full detail) would not provide a large gain, since, anyway local approximations cannot lead to exact results. For the same reasons we also neglect inertial corrections to the boundary conditions in (6) and in the chemical potential (7).

However, we have checked that a model with a continuous transition in azimuthal stress (replacement of $`\mathrm{sin}^4\theta `$ by $`\mathrm{sin}^2\theta `$ in Eq. 9) gives qualitatively the same results.

Now Eq. (3), together with (1) and (9), is a closed third order differential equation for the shape $`r(\theta )`$ which can be easily integrated numerically. It provides both the crack shape (Fig. 1) and a selection of $`v/v_R`$ and $`r_0`$ as functions of the dimensionless driving force $`\mathrm{\Delta }=K^2(1\nu ^2)/2E\alpha `$. The results are given in Figures 2 and 3.

One of the main results is that the upper limit for the steady state crack velocity is appreciable below the Rayleigh speed, as known from experimental results. The instantaneous velocity in the non steady state regime can of course reach higher values Sharon99 ; Fineberg99 . For relatively low driving forces, the growth velocity increases with increasing $`\mathrm{\Delta }`$, but for higher values of $`\mathrm{\Delta }`$, it even decreases. Simultaneously, the tip becomes sharper at first, but then blunts again.

Stability — Although the decrease of the velocity $`v`$ as function of $`\mathrm{\Delta }`$ might be naively understood as a sign of instability, the model itself is stable: We performed a straightforward but tedious numerical stability analysis and found no unstable modes. The point is, that only $`v/v_R`$ decreases, but $`vr_0^3/D`$ increases. Only a decrease of this parameter, which appears in combination with the dissipative coefficient $`D`$, would be a real sign of instability.

Nevertheless, our solution is subject to the ATG instability above a critical threshold of the driving force $`\mathrm{\Delta }`$. The de-stabilizing effects stem from the non-local elastic contributions which we neglected in our model, but which are obviously present in the real problem. This can be proved by simple dimensional analysis arguments: The characteristic wavelength of the ATG instability is $`\lambda E\alpha /\sigma _{\tau \tau }^2`$ Asaro72 ; in the tip region it reads $`\lambda ^{(tip)}r_0/\mathrm{\Delta }`$. Thus, as soon as a certain critical driving force $`\mathrm{\Delta }_c`$ is exceeded, the characteristic wavelength of instability fits into the tip region. The material independent number $`\mathrm{\Delta }_c`$ is the threshold for the instability of the steady state solution, whereas $`\mathrm{\Delta }=1`$ is the Griffith point. Since, according to the steady state solution, $`v/v_R`$ is a universal function of the dimensionless parameter $`\mathrm{\Delta }`$, the threshold of instability in terms of $`v/v_R`$ is also essentially material independent. It is important to note that our steady state predictions are valid only below the threshold of instability. Thus, the main part or all of the decrease of velocity versus $`\mathrm{\Delta }`$ is screened by the instability. Beyond the instability point, the behavior of the system is governed by the full time dependent evolution. Hence we expect the ATG instability to be the relevant mechanism for the experimentally observed microbranching instability Sharon99 ; Fineberg99 . In contrast to the long wave instability Brener98 , this instability is localized in the tip region and cannot be suppressed by convective effects.

Conclusion — We have developed a self-consistent continuum model for crack propagation in homogeneous media. Both ingredients of our theory, the linear elasticity which is valid everywhere in the bulk and surface diffusion which provides a mass transport and dissipative mechanism for crack propagation, are well established. The model is essentially parameter free, leading to the prediction that two dimensionless quantities, the crack velocity - $`v/v_R`$ and crack tip radius - $`r_0(v_R/D)^{1/3}`$ are universal functions of the dimensionless driving force $`\mathrm{\Delta }`$. Strictly speaking, these functions still depend weakly on the Poisson ratio. We note that these statements, together with the prediction of the tip instability above some critical velocity, are based on the general structure of our theory and do not involve the specific modeling of the surface stresses. The specific results given in the figures should differ from exact solutions only quantitatively.

It is important to realize that our model does not contradict classical theories Freund90 , but contains more information. This allows to calculate both the crack velocity and the fracture energy, while the classical theories predict only a relation between these quantities (the integral energy balance). In our model, this energy balance is not fulfilled by approaching the Rayleigh speed limit with increasing driving force, as the classical theories predict, but by an increase of the fracture energy by tip blunting, which eventually leads to the tip instability.

The scale of velocity is set by $`v_R`$ and is independent on the diffusion coefficient. However, the length scale of the tip is set by $`(D/v_R)^{1/3}`$. As we already noted, strong dissipation brings the local crack tip temperature close to the melting temperature. The diffusion coefficient $`D_s`$ is then about $`10^7\mathrm{m}^2\mathrm{s}^1`$ Neumann72  and independent of outside temperature. Consequently, the lengthscale of the tip is of order of atomic units. We hope, that even for these small scales our general, qualitative predictions remain correct.

###### Acknowledgements.

We thank S. Iordanskii, V. I. Marchenko, H. Müller-Krumbhaar, H. Schober, D. Temkin and H. Trinkaus for useful discussions. This work has been supported by the Deutsche Forschungsgemeinschaft (Grant SPP 1120).