text:
# Molecular Signatures from Gene Expression Data

## 1 Introduction

“Molecular signatures” or “gene-expression signatures” are a key feature in many studies that use microarray data in cancer research (e.g., Alizadeh et al.,, 2000; Golub et al.,, 1999; Rosenwald et al.,, 2002; Shipp et al.,, 2002; Pomeroy et al.,, 2002). Shaffer et al., (2001, p. 375) refer to signatures as “(…) genes that are coordinately expressed in samples related by some identifiable criterion such as cell type, differentiation state, or signaling response” (emphasis is ours). Molecular signatures are often used to model patients’ clinically relevant information (e.g., prognosis, survival time, etc) as a function of the gene expression data, but instead of using individual genes as predictors, the predictors are the signature components or “metagenes”.

If we are successful searching for a signature, then we will be able to model, for instance, the probability of developing a metastasis as a function of a few signature components or metagenes where each signature component is made of genes that show strong coexpression. Thus, molecular or gene expression signatures can be important both for diagnostic purposes and for providing information about the biological mechanisms underlying certain conditions by highlighting genes that both coexpress and are related to that condition.

In spite of the widespread use of the term “molecular signature”, no explicit definition is available. Following the conventions of the literature (e.g., Alizadeh et al.,, 2000; Golub et al.,, 1999; Rosenwald et al.,, 2002; Shipp et al.,, 2002; Pomeroy et al.,, 2002; West et al.,, 2001), and building upon the definition above (Shaffer et al.,, 2001, p. 375), we will consider a signature to be composed of one or more signature components or metagenes, where each signature component is a weighted combination of one or more coexpressed genes, and such that statistical models that use signatures both have good predictive performance and are easy to interpret biologically. Interpretation is made easier because the prediction is based on signature components that are weighted averages of subsets of tightly coexpressed genes, which can help when attempting to relate specific biological features to, for example, particular alterations on a metabolic pathway. Based upon the above references, we can try to formalize these goals by requiring that signatures and signature components satisfy the conditions shown in Figure 1.

\[Figure 1 about here\]

The conditions in Figure 1 reflect a very specific biological model. Our objective is to develop a statistical method appropriate for this biological model. By using a method that tries to fulfill those conditions we can also provide evidence that, for any particular case, the underlying biological assumptions behind this attempt are inconsistent with the data or, in other words, the the assumptions embodied in Figure 1 are inappropriate. As will be discussed later, our method is an attempt to map a particular biological model into a statistical method, but other statistical approaches would be more appropriate if more complex biological models are regarded as appropriate.

### 1.1 Limitations of alternative methods

A variety of approaches have been used to identify molecular signatures. A review is provided in the supplementary material. Briefly, most methods return either a single signature component (e.g., Golub et al.,, 1999; Hedenfalk et al.,, 2001; van Belle,, 2002) which is a weighted average of a set of genes, or several signature components (e.g., West et al.,, 2001; Ramaswamy et al.,, 2001; Hastie et al., 2001a, ; Huang et al., 2003b, ; Antoniadis et al.,, 2003) which are often obtained using dimension reduction techniques (e.g., principal component analysis \[PCA\], partial least squares \[PLS\], sufficient dimension reduction) either on the complete set of genes or on a preselected subset.

The most common problems of available methods are:

* Genes within signature components do not necessarily show tight coexpression: no method makes tight coexpression a requirement to be fulfilled.
* The interpretation of components is very difficult for most methods that use PCA or PLS, since all the genes to which PCA or PLS is applied have loadings on each component.
* The search for components in many PCA or clustering of genes methods is carried out without incorporating information from the dependent variable.
* Most methods are designed for a specific type of task (e.g., classification or survival, but not both) and would be difficult to extend to other types of dependent variables.

Our objective in this paper is to propose a method that overcomes these problems. Specifically, our method returns signature components of tight coexpression (and thus, signature components that should ease interpretation) with good predictive performance. Although in this paper we focus on class prediction, our method can be used with different types of dependent data (continuous, categorical, survival), and thus sets up a general framework for finding gene expression signatures regardless of the type of dependent variable.

Based on the proposed operational definition of signatures (see figure 1), we first discuss the key elements of our proposed method. Next, we evaluate its predictive performance and finally with discuss it relation with other methods and problems of biological interpretability. Further details about the algorithm, evaluation of recovery of signatures, and a longer review of alternative approaches are provided in the supplementary material.

## 2 Methods

### 2.1 Key elements of the proposed method

Our objective is to directly fulfill the conditions in figure 1. We start our search with a seed gene that will be the skeleton of a signature component; this first signature component is found so that genes within the component show tight coexpression and the prediction error is acceptable. We repeat this process (find seed gene for a component and then obtain the whole component) greedily, until no further components are needed. The main steps of the algorithm are shown in Figure 2. In this section we explain how the conditions in figure 1 can be fulfilled and provide a geometrical interpretation of the algorithm.

\[Figure 2 about here\]

#### 2.1.1 Fulfilling signature requirements

A common and simple way of characterizing a signature component is to use linear combinations (weighted averages) of the genes that belong to that signature component (e.g., West et al.,, 2001; Shipp et al.,, 2002; Ramaswamy et al.,, 2003; Rosenwald et al.,, 2002). Although we could characterize a signature component using several different linear combinations of the genes of that component, most methods (but see Liu et al.,, 2002) characterize a signature component or metagene using only one linear combination. A single metagene per signature component simplifies interpretation, and is implicit in the requirement that each gene of a signature components should show a strong correlation with the signature component.

Thus, to fulfill requirement one in figure 1, we can use Principal Component Analysis (PCA —which is closely related to Singular Value Decomposition \[SVD\]). PCA yields “the best” representation (or “least distorting” representation, in the least squares sense) of the original data (e.g., Jolliffe,, 2002; Krzanowski,, 1998; Morrison,, 1990) in a subspace of reduced dimensions. The first PC is the best 1-dimensional representation of the original genes of this signature component. If the genes of the signature component are tightly coexpressed, then each of these genes should show a high correlation with the signature component, as we required above (this will also mean that the percentage of variance in the original gene expression data explained by this first PC will be high —see supplementary material). After running the procedure, each signature component will be made of tightly coexpressed genes (we require that all the genes in a component show a correlation larger than a pre-specified threshold of $`r_{min}`$).

In contrast to some previous methods which use PCA or PLS over the complete set of genes, there is no need for our method to return components which are uncorrelated or orthogonal: there is no biological argument that requires that signature components be orthogonal, uncorrelated, or independent (see discussion). For ease of interpretation we will additionally require that no gene belongs to more than one signature component. (In other words, each gene in the original data matrix belongs to either one, and only one, signature component or to none.)

The second and third requirements of figure 1 can be incorporated by adding new signature components only if they result in a relevant reduction of prediction error, and retaining genes in a signature unless they produce large increases in prediction error. In other words, we will penalize adding signature components, but will try to obtain large signature components. The reason is that the search for molecular signatures is often pursued to provide biological insights into coexpressed genes related to conditions and, therefore, minimization of prediction error is not the only goal. Thus, if there are potential trade-offs between prediction error and “biologically interpretable signatures”, the researcher should have the option of modifying the terms of this trade-off flexibly.

#### 2.1.2 Searching for signature components and a geometrical interpretation

Our objective is, thus, to maximize predictive performance using signature components that satisfy that the correlation of each gene in a signature component with the signature component is larger than a given threshold. However, the discussion so far does not indicate how to find the signature components and, given the dimensionality of the problem, an exhaustive search for the optimal solution is not possible. Since we require that each component be highly correlated with the genes of that component, we can start the search with genes that have good predictive abilities on their own. Once we find an initial “seed gene”, we build an initial candidate signature component by including all “promising genes” (e.g., all those with a minimum correlation with the seed gene), and later reduce the signature component eliminating genes until the conditions of minimum correlation with 1st PC (all genes have a correlation with the 1st PC $`>r_{min}`$) and predictive performance are met. (If this elimination eliminates all genes except the seed gene, then, of course, the two requirements are met).

The method proposed here is a heuristic search that has an intuitive geometrical interpretation. We require that each component be highly correlated ($`>r_{min}`$) with the genes in the component, which is equivalent to saying that the vector of the component must have a similar direction as the vectors of each gene in variable space (the space where subjects are the axes). Therefore, no matter which genes belong to a signature component, the component will have a similar direction as any of its genes. Then, it seems reasonable to start the search with the direction that has the best predictive ability, the seed gene; this seed gene is the single direction in space that most contributes to separation of the groups in a classification problem; analogous for regression or survival analysis. When we form the complete signature component, all other genes of the signature component have directions that are similar to that of the seed gene. Together, all the genes of a signature component move the direction slightly, but this shift is possibly towards directions that contribute more to separation of groups (or that at least do not degrade the separation too much) and never moves us far away from the original seed gene. This process is repeated until the addition of new signature components does not achieve any relevant decrease in prediction rate, or until a maximum pre-specified number of signature components is reached. The algorithm is shown in figure 2. Further details are given in the supplementary material.

#### 2.1.3 Choice of underlying classifier

In this paper we will be dealing with a classification problem. Each signature component is used as a predictor variable for a classifier. Of the available classification methods, we have used DLDA (diagonal linear discriminant analysis), a version of linear discriminant analysis which assumes the same diagonal variance-covariance matrix for all the classes (Dudoit et al.,, 2002), and NN (k-nearest neighbor, with k = 1), a simple non-parametric rule that assigns a test sample to the class of the closest training sample (where closeness is measured using Euclidean distance in the space whose dimensions are the signature components). KNN and DLDA have been repeatedly shown to perform as well as, or better than, many competing methods with microarray data (Dudoit et al.,, 2002; Romualdi et al.,, 2003). In addition, DLDA and KNN are simple to implement and interpret. Dudoit et al., (2002) used an adaptive procedure to estimate the optimal number of neighbors to use with KNN; that can be time consuming, and we have fixed K = 1, since this is often a successful rule (Ripley,, 1996; Hastie et al., 2001b, ). As discussed in the supplementary material, other classifiers can be used.

## 3 Comparing predictive performance with established methods

Here we compare the predictive performance of our method with that of three well established methods, support vector machines, KNN, and DLDA, using several “real data” sets. The supplementary material reports simulation studies that show that the suggested method can indeed recover signatures when they are present in the data.

Predictive performance is evaluated using 10-fold cross-validation (i.e., the complete algorithm shown in Figure 2 is applied to each of the 10 “training sets”). This 10-fold cross-validation was repeated 20 times under each condition. The error rates shown are not the CV error rates obtained in steps 1 and 5 of the signature algorithm (see Figure 2), since those are biased down (see supplementary material); the error rates shown are the error rates obtained from cross-validating the complete procedure.

Because an important parameter of our method might be $`r_{min}`$, the minimal absolute correlation between each gene in a signature component and the signature component, we have evaluated the performance of the signature method using a set of values of $`r_{min}`$ that covers a “biologically interesting” range: $`\{.60,.65,.70,.75,.80,.85,.90,.95\}`$. In addition, we have also examined the differences between using $`c_1=c_2=1`$ compared to $`c_1=c_2=0`$; the first corresponds to the usual “1 se rule” and should lead to more interpretable results ($`c_1`$ and $`c_2`$ are related to how much we penalize adding a new signature component and how much we penalize eliminating genes from signature components; see details in supplementary material).

### 3.1 The data sets

From Golub et al., (1999). The original data, from an Affymetrix chip, comprises 6817 genes, but after filtering as done by the authors we are left with 3051 genes. Filtering and preprocessing is described in the original paper and in Dudoit et al., (2002). We used the training data set of 38 cases (27 ALL and 11 AML) in the original paper (the observations in the “test set” are from a different lab and were collected at different times). This data set is available from \[http://www-genome.wi.mit.edu/cgi-bin/cancer/datasets.cgi\] and also from the Bioconductor package multtest (\[http://www.bioconductor.org\]).

From Ramaswamy et al., (2003). We used the data from the 12 metastatic tumors and 64 primary tumors. The original data set included 16063 genes from Affymetrix chips. The data (DatasetA\_Tum\_vsMet.res), downloaded from \[http://www-genome.wi.mit.edu/cgi-bin/cancer/\], had already been rescaled by the authors. We took the subset of 9376 genes according to the UniGene mapping, thresholded the data, and filtered by variation as explained by the authors. The final data set contains 9868 clones (several genes were represented by more than one clone); of these, 196 had constant values over all individuals.

From Ross et al., (2000). The data, from cDNA arrays, can be obtained from \[http://genome-www.stanford.edu/sutech/download/nci60/index.html\]. The raw data we used, which is the same as the data used in Dudoit et al., (2002); Dettling & Bühlmann, (2003), is the one in the file “figure3.cdt”. As in Dudoit et al., (2002); Dettling & Bühlmann, (2003) we filtered out genes with more than two missing observations and we also eliminated, because of small sample size, the two prostate cell line observations and the unknown observation. After filtering, we were left with a 61 x 5244 matrix, corresponding to eight different tumor types (note that, as done by previous authors, we did not average the two observations with triplicate hybridizations). As in Dudoit et al., (2002) we used 5-nearest neighbor imputation of missing data using the program GEPAS (Herrero et al.,, 2003) (http://gepas.bioinfo.cnio.es/cgi-bin/preprocess); unlike Dudoit et al., (2002), however, we measured gene similarity using Euclidean distance from the genes with complete data, instead of correlation: Troyanskaya et al., (2001) found Euclidean distance to be an appropriate metric. Finally, as in (Dudoit et al.,, 2002, p. 82) gene expression data were standardized so that arrays had mean 0 and variance 1 across variables (genes).

From van Belle, (2002). The data, from Affymetrix arrays, were downloaded from \[http://www.rii.com/publications/2002/vantveer.htm\] (we used the files ArrayData\_less\_than\_5yr.zip, ArrayData\_greater\_than\_5yr.zip, ArrayData\_BRCA1.zip, corresponding to 34 patients that developed distant metastases within 5 years, 44 that remained disease-free for over 5 years, and 18 with BRCA1 germline mutations and 2 with BRCA2 mutations). As did by the authors, we selected only the genes that were “significantly regulated” (see their definition in the paper and supplementary material), which resulted in a total of 4869 clones. Because of the small sample size, we excluded the 2 patients with the BRCA2 mutation. We used 5-nearest neighbor imputation for the missing data, as for the NCI 60 data set. Finally, we excluded from the analyses the 10th subject from the set that developed metastases in less than 5 years (sample 54, IRI000045837, in the original data files), because it had 10896 missing values out of the original 24481 clones, and was an outstanding outlying point both before and after imputation. The breast cancer dataset was used both for two class comparison (those that developed metastases within 5 years vs. those that remain metastases free after 5 years) and for three group comparisons.

Therefore, we use three datasets in which the problem is classification into two classes (leukemia, adenocarcinoma, breast cancer), one dataset with a three class problem (breast cancer) and one dataset with an eight class problem.

### 3.2 The competing methods

We have used three methods that have shown good performance in reviews of classification methods with microarray data (Dudoit et al.,, 2002; Romualdi et al.,, 2003).

DLDA is the maximum likelihood discriminant rule, for multivariate normal class densities, when the class densities have the same diagonal variance-covariance matrix (i.e., variables are uncorrelated, and for each variable, its variance is the same in all classes). This yields a simple linear rule, where a sample is assigned to the class $`k`$ which minimizes $`\mathrm{\Sigma }_{j=1}^p(x_j\overline{x}_{kj})^2/\widehat{\sigma }_j^2`$, where $`p`$ is the number of variables, $`x_j`$ is the value on variable (gene) $`j`$ of the test sample, $`\overline{x}_{kj}`$ is the sample mean of class $`k`$ and variable (gene) $`j`$, and $`\widehat{\sigma }_j^2`$ is the (pooled) estimate of the variance of gene $`j`$ (Dudoit et al.,, 2002). In spite of its simplicity and its somewhat unrealistic assumptions (independent multivariate normal class densities), this method has been found to work very well.

KNN is a non-parametric classification method that predicts the sample of a test case as the majority vote among the k nearest neighbors of the test case (Ripley,, 1996; Hastie et al., 2001b, ). To decide on “nearest” here we use, as in Dudoit et al., (2002), the Euclidean distance. The number of neighbors used (k) is chosen by cross-validation as in Dudoit et al., (2002): for a given training set, the performance of the KNN for values of $`k`$ in $`\{1,3,5,\mathrm{},21\}`$ is determined by cross-validation, and the $`k`$ that produces the smallest error is used.

SVM are becoming increasingly popular classifiers in many areas, including microarrays (Furey et al.,, 2000; Lee & Lee,, 2003; Ramaswamy et al.,, 2001). SVM (with linear kernel, as used here) try to find an optimal separating hyperplane between the classes. When the classes are linearly separable, the hyperplane is located so that it has maximal margin (i.e., so that there is maximal distance between the hyperplane and the nearest point of any of the classes) which should lead to better performance on data not yet seen by the SVM. When the data are not separable, there is no separating hyperplane; in this case, we still try to maximize the margin but allow some classification errors subject to the constraint that the total error (distance from the hyperplane in the “wrong side”) is less than a constant. For problems involving more than two classes there are several possible approaches; the one used here is the “one-against-one” approach, as implemented in libsvm Chang & Lin, (2003). Reviews and introductions to SVM can be found in Burgues, (1998); Hastie et al., 2001b .

For each of these three methods we need to decide which of the genes will be used to build the predictor. Based on the results of Dudoit et al., (2002) we have used the 200 genes with the largest $`F`$-ratio of between to within groups sums of squares. Dudoit et al., (2002) found that, for the methods they considered, 200 genes as predictors tended to perform as well as, or better than, smaller numbers (30, 40, 50 depending on data set).

We evaluated predictive performance using 10-fold cross-validation; the results shown are from 20 replications of the 10-fold cv process. In all cases, cross-validation includes gene selection (Ambroise & McLachlan,, 2002; Simon et al.,, 2003); in other words, for the three competing methods and the signature algorithm the selection of genes is carried out within each of the 10 “training sets” of the cross-validation. Thus, we insure that the subjects for which prediction is performed have not been used for the gene selection process.

## 4 Results

Predictive performance is shown in Figure 3. Predictive performance changes very little from using DLDA vs. NN, or setting $`c_1=c_2=1`$ vs. $`c_1=c_2=0`$; figures for four combinations of classifier and values of $`c_1,c_2`$ are shown in the supplementary material. Comparing $`c_1=c_2=1`$ with $`c_1=c_2=0`$ (see Figures in supplementary material) does not show any relevant differences in predictive performance; of course, there are differences in the outcome because, not surprisingly, using $`c_1=c_2=0`$ tends to result in more signature components of smaller numbers of genes per component, and higher correlations between components.

\[Figure 3 about here\]

\[Table 1 about here\]

With respect to $`r_{min}`$, except for the NCI data set, and slightly for the Breast cancer with 3 classes data set, changes in $`r_{min}`$ have little effects on predictive performance. In Table 1 we show the median number of components, median total number of genes in a signature, and median average number of genes per component obtained in 200 bootstrap runs using $`c_1=c_2=1`$ with $`r_{min}=0.85`$ and $`r_{min}=0.6`$, and using NN as the classifier (see also section 4.1). It can be seen that changes in $`r_{min}`$ do affect the outcome in terms of number of number of genes per component. These results seem to indicate that choice of $`r_{min}`$ can probably be guided more by interpretability concerns (whether we want larger signature components of looser coexpression or smaller signature components of tight coexpression) than by concerns over predictive ability.

Finally, the performance of the signature method is only slightly worse than that of the three competing classifiers, except for the NCI data set. As seen in Table 1, most of the signatures, specially with $`r_{min}=0.85`$, used very few components of very few genes each (compared to the use of 200 genes for the competing classifiers); thus, the predictive performance is achieved using a very small number of genes and thus potentially facilitating a simple biological interpretation. In the case of the NCI data set, there are eight classes with only 61 samples, and in most cases the signature component only returned between 1 and 3 components. Probably the forward sequential addition of components in the signature method has affected negatively the predictive capabilities, because any single addition was most likely incapable of resulting in a large enough decrease of prediction error to justify further addition of components (in contrast to using, directly, 200 genes as predictors).

### 4.1 Stability of results

To evaluate the stability of results, we rerun the complete procedure on all data sets using the bootstrap (Efron & Tibshirani,, 1993; Davison & Hinkley,, 1997) with 200 bootstrap samples, similar to what Efron & Gong, (1983) do to evaluate a complex fitting procedure. We run the procedure for settings of $`c_1=c_2=1`$ with $`r_{min}=0.85`$ and $`r_{min}=0.6`$, and using NN as the classifier. The results are shown in Table 2. The bootstrap results indicate that, when using a high value of $`r_{min}`$ we rarely obtain similar solutions repeatedly; some data sets, however, seem to yield more stable solutions (e.g., Leukemia data sets) and, not surprisingly, if the $`r_{min}`$ criterion is set to less stringent values, results tend to be more repeatable.

\[Table 2 about here\]

## 5 Discussion

### 5.1 Similarities and differences with other methods

Our method is unique because it simultaneously searches for sets of genes that are tightly coexpressed and will lead to good predictive performance. The search for the sets of genes is carried out using the information from the dependent variable (at the first stage —when selecting the seed gene—, and at the pruning stage of reducing the signature component —when genes that lead to decreased predictive performance are eliminated from the signature component —; see figure 2).

One important difference between our proposed method and most previous approaches that use PCA is that, by performing PCA only on subsets of genes, our method returns signature components where genes show tight coexpression. Because returned components are not orthogonal, and simple components are an explicit goal, our approach is actually closer to some ideas implemented in SAS’s PROC VARCLUS, which is similar to factor analysis with oblique rotation and can be used to obtain clusters of variables, of relatively simple interpretation, to be further used in model building SAS Insitute, (1999); Nelson, (2001); see also Harrell, (2001).

Using PCA on subsets of genes, instead of the complete set of genes is crucial because it makes interpretation easier and allows for subsets of tight coexpression. Simple PCA and related methods (Rousson & Gasser,, 2003; Vines,, 2000; Jolliffe,, 2002) as well as SAS PROC VARCLUS also try to achieve components of tight coexpression, but many of these approaches cannot be applied with $`pn`$, and all of them carry out the PCA without using the information from the dependent variable, which in our case is a fundamental requirement since the sets of genes with tightest coexpression would be irrelevant for our purposes if they are not related to the dependent variable we are trying to model. This difference in objectives is also evident because our aim is not to explain the most variance in the genes (as in most simple PCA approaches) nor maximize variable explained across all clusters (such as in PROC VARCLUS). Overall summarization of information is not important for our problem, because we are interested in prediction, and we often have the suspicion that most of the genes in the array are not related to the outcome variable. Finally, our approach can result in signature components which are correlated (sometimes strongly), but this is not inherently a problem because there is no biological reason to suggest that the underlying biological causes or factors ought to be independent or uncorrelated; moreover, if the true underlying causes are not orthogonal, using a method such as PCA can lead to interpretational and conceptual difficulties because each biological cause will be spread over several orthogonal components (Houle et al.,, 2002): non-orthogonal biological causes are inconsistent with procedures such as PCA and PLS.

Bayesian classification trees using the 1st PC from gene clusters (Huang et al., 2003a, ) and block PCA (Liu et al.,, 2002) also used PCA on subsets of genes, instead of the complete set of genes. In both cases their subsets of genes were obtained using criteria that did not make any use of the information from the dependent variable. In addition, in Huang et al., 2003a  the metagenes are not necessarily of subsets of tightly coexpressed genes. In Liu et al., (2002) there is an explicit criterion of % variance accounted for, but often the number of components used to summarize a subset of genes is too large to allow for easy interpretation (11 to 16 principal components).

Supervised harvesting of expression genes (Hastie et al., 2001a, ) also works with clusters of subsets of genes, which are then used in a predictive model. As before, however, the clustering is carried out without using information from the dependent variable; even if the selection of which subsets or clusters to use in the model uses the information from the dependent variable, the very first step of clustering genes does not, and can therefore be unable to recover sets of tightly coexpressed genes that are good predictors. Finally, the “Wilma” and “Pelora” methods (Dettling & Bühlmann,, 2002, 2004) do use the information from the dependent variable in the formation of clusters of genes; however, there is no explicit objective of achieving tight gene coexpression within clusters and thus, how tightly coexpressed genes are in each metagene cannot be specified in advance; in addition, Wilma weights each gene equally within a cluster (only possible weights are +1 and -1) whereas we use PCA on unscaled genes, thus allowing genes to play a different role in the specification of the direction of the signature component (genes with larger among-subject variance play a more important role in determining direction).

The rest of the alternative methods differ strongly from our proposed method, either because they do not return subsets of genes but components with loadings from all genes (e.g., PLS based methods), or return subsets of genes where there is not requirement of tight coexpression (e.g., weighted gene voting, Golub et al.,, 1999).

After gene selection and dimension reduction (i.e., the use of only the 1st PC, that collapses all the information from the genes of a signature component onto one dimension), the predictive model of our choice is fitted. In this sense, the method as presented here is “just” a DLDA or NN that uses signature components instead of genes as the predictors. The choice of DLDA and NN was made based on published results that showed their excellent performance with microarray data. In particular, Dudoit et al., (2002) showed that other forms of discriminant analysis tended to perform much worse because of the small ratio sample size/parameters needed to estimate covariances and different variances per group; however, since in most cases our method returns just a few signature components, other types of discriminant analysis that use the information from the covariance of the predictors (e.g., linear discriminant analysis and quadratic discriminant analysis) might prove useful.

### 5.2 Coexpression: across-group and within-group

The algorithm can include in a signature component genes that show no correlation within groups but that show correlation among groups because they are far apart in the multidimensional space, and the correlation coefficient is computed across the whole, pooled, sample. The algorithm might even include in the same signature component genes that have very different patterns of correlation in different groups, if they still show sufficiently strong correlation over the pooled sample.

To our knowledge, this issue has not been explicitly addressed in any other approach to the signature problem (see reviews in supplementary material). However, probably the most biologically relevant components are those where there is strong correlation within groups, because this is a more reliable indicator of coordinated expression.

A possible solution is, for example, to only accept results for a signature component if a principal component analysis over the pooled sample after centering the data with respect to the group means yields a relevant first eigenvalue; for added robustness, we might want to use the trimmed mean. This approach, however, does not directly address if there are different multivariate orientations in different groups of subjects, and how these orientations within-group relate to the across-group orientation. In particular, the case where several groups not only have the same first principal component, but are lined up along a common axis, known as “allometric extension” (Hills,, 1982; Bartoletti et al.,, 1999), might constitute the most natural type of signature component.

Krzanowski (see reviews and summary in Krzanowski,, 1998; Jolliffe,, 2002) has proposed a method to directly compare the subspaces defined by the principal components of each of the groups. In our case, as we only use the first principal component of a set of genes to define a signature component, for each signature component we can compare the first principal component of each group. An example using the NCI 60 data set is shown in the supplementary material. This method only compares the orientation of the principal components (the eigenvectors) but does not compare the location of the multivariate means. The EDDA (Bensmail & Celeux,, 1996) and common principal components (Flury,, 1988) approaches provide frameworks to examine differences among covariance matrices of particular interest in the context of discrimination among groups. Nevertheless, it must be remembered that biological interpretation of results is not always straightforward, specially when the underlying biological factors are not orthogonal (Houle et al.,, 2002), and that we need to asses the need for robust methods (Boente et al.,, 2002) and power issues related to the small sample sizes common in microarray data. We are currently investigating some of these issues; however, in the presence of the difficulties associated with interpretation, low statistical power, and possible lack of robustness to outlying observations, right now the best recommendation could be to conduct an examination of patterns of within vs. among coexpression for the returned signatures.

### 5.3 Stability and biological relevance of signatures

Our method follows from an operationalization of signatures and signature components (figure 1). If the ideas embodied in figure 1 have empirical support, then it might be possible to build good predictive models using just a few, very specific biological features that can be related to alterations on particular pathway. Our method is unique in this regard because it returns only a few signature components where genes within each component are tightly coexpressed.

On the other hand, our method can indicate that the data are inconsistent with the ideas behind figure 1: if the signature components show high instability, this is evidence that very different models can be obtained from the data. Widely different models from a reasonably large data set cast doubt on the idea that a few, easily interpretable, signature components strongly correlated with the expression of a few key genes are associated to the clinical outcome of interest. In the context of building predictors from gene expression data, Somorjai et al., (2003) (see also Zhang et al.,, 2003) have emphasized that this non-uniqueness leads to interpretational difficulties and should make researchers skeptical about the biological relevance of any set of predictors; moreover, they explain how this non-uniqueness can arise from dataset sparsity. Of course, both of these issues are relevant to the present proposed methodology. None of the data sets examined in this paper yield stable signatures when we use stringent criteria of gene coexpression (see 4.1 and Table 2 with $`r_{min}=0.85`$), although some of the data sets are somewhat stable from run to run when the $`r_{min}`$ is set to small values (but other data sets show signatures that vary widely from run to run). These results add to the above references in the sense that biological interpretation should be carried out very cautiously, and emphasize the difference between attempting to build good predictors and attempting interpretation (see also Breiman,, 2001). More relevant to the current work, the present results indicate that simple models of molecular signatures warrant further critical scrutiny, and that it might be extremely hard to identify molecular signatures from such sparse data sets (see Rhodes et al.,, 2004, for a meta-analysis attempt to identify stable “signatures”).

We must recognize that these results are preliminary for two main reasons. First, establishing that two or more signature components are different probably requires additional information besides the identities of the genes; for instance, information from Gene Ontology, or known participation on certain regulatory networks. Second, we have used two different $`r_{min}`$ thresholds, but it is unclear what constitute “biologically reasonable” patterns of covariation between genes that are to belong to the same signature component. Nevertheless, even if preliminary, these results emphasize the need for further work in the operationalization and explicit definition of what we mean by molecular signatures, careful consideration of the stability of results, and critical assessment of the sample sizes need to reliably identify molecular signatures.

### 5.4 Alternative statistical methods for alternative biological models

As mentioned in the introduction, we started by trying to clarify, conceptually, what is often understood by molecular signature (see Figure 1), and then devised a statistical method to fulfill those requirements. The biological model underlying the suggested method is one where most of the genes are not relevant for prediction, relevant genes are involved in one and only one signature component (i.e., non-overlapping signature components), and the signature components are common, and have similar covariance matrices, in different groups.

However, other biological models are plausible, and for those biological models other statistical methods would be more appropriate. The simultaneous clustering and classification approach in Jörnsten & Yu, (2003) could be extended by placing restrictions on the covariance matrix (i.e., require a minimum correlation between genes) but possibly allowing for different covariance matrices among groups; thus, we could address directly issues of different across vs. within-group correlations (see section 5.2), within a formal inferential framework. Biological models where signature components are not common and/or do not behave similarly in different groups could be investigated using modifications of the Plaid model of Lazzeroni & Owen, (2002) (see also Turner et al.,, 2004). In addition, genes with the highest correlation need not be the best candidates for being in the same biological pathway; activity in a pathway might just require that precursor genes get activated, but once a threshold is reached, it might not be very important by how much the threshold is exceeded. This type of behavior could preclude strong correlation between genes that belong to the same pathway. This can be modelled building upon the latent class methods of Parmigiani and colleagues (Parmigiani et al.,, 2002; Garrett & Parmigiani,, 2003; Scharpf et al.,, 2003), where signature components are based on under-, over- or baseline expression (instead of expression levels). Work along these lines is currently in progress in our group.

### 5.5 “Just” dimensional reduction?

Even if the current method fails to identify stable features that can be associated to molecular signatures, it can be a useful dimension reduction tool. Difficulties associated with a simple mapping of the returned “signature components” to pathways, and problems derived from instability of the found components also affect any other of the existing alternative methods. Thus, in the presence of instability of results, it is more appropriate to regard this method as dimension reduction tool that could lead to simple biological interpretation. The simple biological interpretation could be helped not only because of the coexpression of the genes that make a signature component, but also because the dimension reduction performed is quite remarkable compared to other methods (the number of signature components and genes returned is very small for the five data sets examined; see Results) with, at most, only a slight decrease in predictive performance. Moreover, the user can control the relative trade-offs between predictive performance and potential interpretability of results (e.g., coexpression of sets of genes) by changing the $`r_{min}`$ parameter (note that if $`r_{min}=1`$ the method becomes essentially either DLDA or NN with forward addition of genes to the model). This flexible modification of the trade-offs between prediction error and interpretability is of great importance in methods that are largely exploratory and oriented towards providing “biologically interpretable” output; in other words, methods for which minimization of prediction error should not be the only goal.

## 6 Conclusions

The most common methods for finding signatures present several deficiencies that do not allow them to return signatures and signature components that fulfill basic biological requirements. After suggesting an operational definition of signature and signature components, we have developed a method that follows directly from what are often considered as the biologically relevant signature characteristics. The method developed returns signature components of tightly coexpressed genes and thus can facilitate biological interpretation. In this paper we have applied the method to classification problems, but this approach in fact sets up a framework that allows us to find signatures regardless of the type of dependent variable. Extension to use other classifiers is straightforward and it should also be easy to incorporate other types of dependent variables to allow, for example, survival analysis. We have also shown that the predictive performance of our method is comparable to that of state of the art methods. Finally, our method not only could facilitate mapping pathological alterations to a few, tightly coexpressed sets of genes, but can also provide evidence that the underlying biological assumptions behind this attempt are inconsistent with the data. In the five data sets analyzed, our results suggest that identification of molecular signatures is questionable under this simple model. These results emphasize the needed for further work on the operationalization of the biological model and the necessity of critical assessment of the stability of putative signatures.

## 7 Acknowledgements

S. Ramaswamy for answering questions about the data and processing steps for the adenocarcinoma dataset, H. Dai for explanations about the breast cancer dataset, and S. Dudoit for answering some questions about the NCI data set. A. Pérez and M. A. Piris for introducing me to, and emphasizing the importance of, molecular signatures. The Bioinformatics Unit at CNIO, C. Lázaro-Perea, Y. Benjamini and F. Falciani for discussion. D. Casado for help with the code. C. Lázaro-Perea and two anonymous reviewers for detailed and thorough comments on the ms. The author was partially supported by the Ramón y Cajal program of the Spanish MCyT (Ministry of Science and Technology); funding partially provided by project TIC2003-09331-C02-02 of the Spanish MCyT.


format: text

text:
# The New Physics at RHIC. From Transparency to High pt Suppression.

## 1 Introduction

With the Relativistic Heavy Ion Collider, RHIC, at Brookhaven National Laboratory a new chapter in the story of the production and study of nuclear matter at extreme energy density and temperature is being written. In collisions between gold nuclei at 100AGeV+100AGeV at RHIC the total energy in the center of mass is almost 40TeV, the largest so far achieved in nucleus-nucleus collisions under laboratory conditions. This energy is so large that if a sizeable fraction of the initial kinetic energy can be converted into matter production, many thousands of particles can be created in a limited volume. It is expected that a region of matter consisting of free quarks and gluons may be formed in the early stages of the collision, the so called Quark Gluon Plasma (QGP). The salient questions to answer are: how much of the initial (longitudinal) kinetic energy is converted to other degrees of freedom, what is the energy density achieved and are there specific signals that carry evidence for the nature of the formed high density state (i.e. is it partonic or hadronic)? In this talk I will try to answer some of these questions based on the results from the first rounds of experiments with RHIC.

RHIC started beam operations in the summer of year 2000 with a short commissioning run colliding Au nuclei at $`\sqrt{s_{NN}}=130`$ GeV). The first full run at top energy ($`\sqrt{s_{NN}}=200`$ GeV) took place in the fall/winter of 2001/2002. The third RHIC run during the winter/spring of 2003 focussed on d+Au and p+p reactions. Here, I concentrate on results from the BRAHMS detector, one of the four major detectors at RHIC, from the two long production runs. BRAHMS is a two arm magnetic spectrometer with excellent momentum resolution and particle identification capabilities for hadrons. The spectrometers subtend only a small solid angle (a few msr) but they can rotate in the horizontal plane about the collision point enabling the collection of data on hadron production over a wide rapidity range (0-4), a unique feature among the RHIC experiments. Details about the BRAHMS detector system may be found in .

## 2 Transparency at RHIC

The stopping of colliding nuclei is a long standing subject. In reactions at energies in the MeV domain atomic nuclei behave like droplets. In gentle collisions, nuclei can fuse forming an excited compound system at rest in the CM frame. In somewhat more violent collisions the droplets may splash and fragment with approximately isotropic distributions in the CM. Even at energies of several GeV pr. nucleon, f.for example in reactions with Au ions at the AGS, where the energy in the nucleon nucleon center of mass system is $`\sqrt{s_{NN}}=5`$ GeV, the colliding nuclei, as seen from the CM frame, experience significant stopping and a large fraction of the initial kinetic energy is converted into excitations which subsequently decay and release energy for example in a direction transverse to the beam direction.

A useful way to quantify the stopping is by the rapidity loss experienced by the baryons in the colliding nuclei. Rapidity is defined as

$$y=\frac{1}{2}ln(\frac{E+p_z}{Ep_z})=\frac{1}{2}ln(\frac{1+\beta cos\theta }{1\beta cos\theta })$$
(1)
where $`E,pz,\beta `$ and $`\theta `$ denote the total energy, longitudinal momentum, velocity and angle relative to the beam axis, respectively, of a particle. If incoming beam baryons have rapidity, $`y_b`$ relative to the CM (which has $`y=0`$) and average rapidity

$$<y>=_0^{y_b}y\frac{dN}{dy}𝑑y$$
(2)
after the collision, the rapidity loss is $`\delta y=y_b<y>`$. Here dN/dy denotes the number of particles per unit of rapidity. Thus, for the extreme case of full stopping: $`\delta y=y_b`$. This corresponds to the situation found at very low energies where all the beam baryons loose all their kinetic energy. In the expression above a complication arises at CM energies large enough to allow for the formation of baryon-antibaryon pairs. Thus the baryon dN/dy distribution to be used is that for the net number of baryons ( i.e. the difference between the number of baryons and antibaryons).

At AGS the number of produced antibaryons is quite small and the the net-baryon distribution is similar to the proton distribution. The net-proton rapidity distribution is centered around $`y=0`$ and is rather narrow. The rapidity loss is about 1 for a beam rapidity of approx. 1.6. At CERN-SPS energies ($`\sqrt{s_{NN}}=17`$ GeV, 158 AGeV Pb+ Pb reactions) the rapidity loss is about 2 for a beam rapidity of 2.9, about the same relative rapidity loss as at AGS. The fact that the rapidity loss is large on an absolute scale means, however, that there is still a sizeable energy loss of the colliding nuclei. This energy is available for particle production and other excitations. Indeed, in collisions at SPS, multiplicities of charged hadrons are about dN/dy=180 around y=0. At SPS another feature is visible (see fig. 1): the net proton rapidity distribution shows a double ’hump’ with a dip around y=0. This is a consequence of two effects: the finite rapidity loss of the colliding nuclei and the finite width of each of the humps, which reflect the rapidity distributions of the protons in the colliding nuclei after the collisions. This picture suggests that the reaction at SPS is beginning to be transparent in the sense that fewer of the original baryons are found at midrapidity after the collisions, in contrast to the situation at lower energies.

We have measured the net proton rapidity distribution at RHIC in the interval $`y=03`$ in the first run with (0-10%) central Au+Au collisions at full energy. The beam rapidity at RHIC is about 5.4. Details of the analysis may be found in refs. and in the recent Ph.D. thesis of P. Christiansen. The results are displayed in fig. 1 together with the previously discussed net-proton distributions measured at AGS and SPS. It is notable that the RHIC distribution is both qualitatively and quantitatively different from those at lower energies. The net number of protons per unit of rapidity around y=0 is only about 7 and the distribution is flat over at least the first unit of rapidity. The distribution increases in the rapidity range $`y=23`$ to an average $`dN/dy12`$. We have not yet completed the measurements at the most forward angles (highest rapidity) allowed by the geometrical setup of the experiment, but we can exploit that there must be baryon conservation in the reactions to try to set limits on the relative rapidity loss at RHIC. This is illustrated in fig. 2, which shows various possible distributions whose integral areas correspond to the number of baryons present in the overlap between the colliding nuclei. From such distributions one may deduce a set of upper and lower limits for the rapidity loss at RHIC. In practice the situation is complicated by the fact that not all baryons are measured. We measure in BRAHMS the direct protons, but only some of the decay protons from for example $`\mathrm{\Lambda }`$. The limits shown in the figure include some reasonable estimates of these effects .

The conclusion is that the absolute rapidity loss at RHIC $`(\delta y=2.2\pm 0.4)`$ is not appreciably larger than at SPS. In fact the relative rapidity loss is significantly reduced as compared to an extrapolation of the low energy systematics . It should be noted that the rapidity loss is still significant and that, since the overall beam energy (rapidity) is larger at RHIC than at SPS, the absolute energy loss increases appreciably from SPS to RHIC thus making available a significantly increased amount of energy for particle creation in RHIC reactions.

## 3 Particle production and energy density

The stopping scenario that we observe at RHIC and which was outlined in the previous section entails that the reaction can be viewed as quite transparent (opaque is perhaps a better word). After the collision, the matter and energy distribution can be conceptually divided up into two main parts, namely a so-called fragmentation region consisting of the excited remnants of the colliding nuclei which have experienced an average rapidity loss of about 2.2 and a central region in which few of the original baryons are present but where significant energy density is collected. This picture is consistent with the schematic one already proposed by Bjorken 20 years ago . The central region (an interval around midrapidity) is decoupled from the fragments. The energy removed from the kinetic energy of the fragments is initially stored in a color field strung between the receeding partons that have interacted. The linear increase of the color potential with distance eventually leads to the production of quark-antiquark pairs. Such pairs may be produced anywhere between the interacting partons leading to an approximately uniform particle production as a function of rapidity. In this picture, the properties of the particle production is also uniform as a function of rapidity (boost invariance). If the density of produced quark-antiquark pairs is sufficiently high, the average distance between them will be low and the binding potential between the colored objects will be small. The objects will become asymptotically free and exist in a plasma like state until the subsequent expansion and lowered density leads to confinement and hadronization.

Figure 3 shows the overall multiplicity of charged particles observed in Au+Au collisions at RHIC  for various collision centralities and as a function of pseudorapidity (pseudorapidity, $`\eta `$, is defined as $`\eta =ln(cot(\theta /2))`$ and is a customary rapidity variable for non identified particles). The figure shows that the multiplicity at RHIC is about $`dN/d\eta =625`$ charged particles pr. unit rapidity around $`\eta =0`$ for central collisions. The production of charged particles in central collisions exceeds the particle production seen in p+p collisions at the same energy by about 40%, when the yield seen in p+p collisions is multiplied by the number of participant nucleon pairs in the overlap region between the colliding nuclei.

Figure 4 shows a recent and more detailed study of the particle production in central collisions as a function of rapidity (the PhD. work of Djamel Ouerdane ). The figure shows the rapidity densities of pions and kaons for central collisions. From such distributions we can construct the ratio of the yields of particles and their antiparticles as a function of rapidity. Figure 5 shows the ratios of yields of antihadrons to hadrons (posititive pions, kaons and protons and their antiparticles). The ratio is seen to be approaching unity in an interval of about 1.5 units of rapidity around midrapidity, suggesting that the particle production is predominantly from pair creation. This is exactly true for pions (ratio of 1), but less so for kaons (ratio= 0.95) and protons (ratio= 0.76). The reason is that there are other processes that break the symmetry between particles and antiparticles that depend on the net-baryon distribution discussed in the previous section. One such process that is relevant for kaons is the associated production mechanism (e.g. $`p+pp+\mathrm{\Lambda }+K^+`$) which leads to an enrichment of positive kaons in regions where there is an excess of baryons. Support for this view is given by fig. 6, which shows the systematics of kaon production relative to pion production as a function of center of mass energy. At AGS, where the net proton density is high at midrapidity, the rapidity density of $`K^+`$ strongly exceeds that of $`K^{}`$. In contrast, at RHIC, production of $`K^+`$ and $`K^{}`$ is almost equal. This situation changes, however, at larger rapidities where the net proton density increases.

Integration of the charged particle pseudorapidity distributions corresponding to central collisions tells us that about 4600 charged particles are produced in each of the 5% most central collisions. Since we only measure charged particles, which are predominantly pions and kaons, as may be seen from fig. 4, and not the neutrals, we multiply this multiplicity by 3/2 to obtain the total particle multiplicity of about 7000 particles.

From the measured spectra of pions, kaons and protons and their antiparticles as a function of transverse momentum we can determine the average transverse momentum for each particle species (fig. 7). This allows us to estimate the initial energy density from Bjorkens formula:

$$ϵ=\frac{1}{\pi R^2\tau }\frac{d<E_t>}{d\eta }$$
(3)
where we can make the substitution $`d<E_t>=<m_t>dN`$ and use quantities from the measured spectral distributions. Since we wish to calculate the energy density in the very early stages of the collision process we may use for R the radius of the overlap disk between the colliding nuclei, thus neglecting transverse expansion. The formation time is more tricky. It is often assumed to be of the order of 1 fm/c, but it might be shorter for heavier or more energetic particles, as the uncertainty relation would tell us. Under these assumptions we find that $`ϵ>5GeV/fm^3`$. This value of the initial energy exceeds by a factor of 30 the energy density of a nucleus, by a factor of 10 the energy density of a baryon and by a factor of 5 the critical energy density for QGP formation that is predicted by lattice QCD calculations.

## 4 Is there thermodynamical and chemical equilibrium at RHIC?

The particle yields measured by BRAHMS also lend themselves to an analysis of the charged particle production in terms of the statistical model ,. Figure 8 shows the ratios of negative kaons to positive kaons as a function of the corresponding ratios of antiprotons to protons for various rapidities at RHIC. The data are for central collisions, and the figure also displays similar ratios for heavy ion collisions at AGS and SPS energies. There is a striking correlation between the RHIC/BRAHMS kaon and proton ratios over 3 units of rapidity. Assuming that we can use statistical arguments based on chemical and thermal equilibrium at the quark level, the ratios can be written

$$\frac{\rho (\overline{p})}{\rho (p)}=exp(\frac{6\mu _{u,d}}{T})$$
(4)
and

$$\frac{\rho (K^{})}{\rho (K^+)}=exp(\frac{2(\mu _{u,d}\mu _s)}{T})=exp(\frac{2\mu _s}{T})\times [\frac{\rho (\overline{p})}{\rho (p)}]^{\frac{1}{3}}$$
(5)
where $`\rho ,\mu `$ and $`T`$ denote number density, chemical potential and temperature, respectively. From equation 4 we find the chemical potential for u and d quarks (often called the baryochemical potential) to be around 25 MeV, the lowest value yet seen in nucleus-nucleus collisions. Equation 5 tells us that for a vanishing strange quark chemical potential we would expect a power law relation between the two ratios with exponent 1/3. The observed correlation is well described by the relationship $`\rho (K^{})/\rho (K^+)=\rho (\overline{p})/\rho (p)^{0.24}`$, i.e. with an exponent that is close to 1/4 suggesting, a finite value of the strange quark chemical potential. A more elaborate analysis for a grand canonical ensemble assuming charge, baryon and strangeness conservation can be carried out by fitting these and many other particle ratios observed at RHIC in order to obtain the chemical potentials and temperatures. An example of such a procedure is shown in fig. 8 and displayed with the full line . Here the temperature is 170MeV. The point to be made here is that the calculation agrees with the data over a wide energy range (from SPS to RHIC) and over a wide range of rapidity at RHIC. This may be an indication that the system is in chemical equilibrium over the considered $`\sqrt{s}`$ and $`y`$ ranges (or at least locally in the various y bins). Separate measurements at RHIC of, for example, elliptical flow also point to local equilibration around midrapidity.

## 5 High pt suppression. The smoking gun of QGP?

The discussion in the previous sections indicates that the conditions for particle production in a interval $`|y|<1.52`$ are radically different than in reactions at lower energies. At RHIC the central zone is nearly baryon free, the considered rapidity interval appears to approximately exhibit the anticipated boost invariant properties, the particle production is large and dominated by pair production and the energy density appears to exceed significantly the one required for QGP formation. The overall scenario is therefore consistent with particle production from the color field, formation of a QGP and subsequent hadronization. But is this interpretation unique and can more mundane explanations based on a purely hadronic scenario be excluded? In spite of the difficulties in reconciling the high initial energy density with hadronic volumes, a comprehensive answer to this question requires the observation of an effect that is directly dependent on the partonic or hadronic nature of the formed high density zone.

Such an effect has recently been discovered at RHIC and is related to the suppression of the high transverse momentum component of hadron spectra in central Au+Au collisions as compared to scaled momentum spectra from p+p collisions . The effect, originally proposed by Bjorken, Gyulassy and others  is based on the expectation of a large energy loss of high momentum partons scattered in the initial stages of the collisions in a medium with a high density of free color charges. According to QCD colored objects may loose energy by radiating gluons by bremsstrahlung. Due to the color charge of the gluons, the energy loss is proportional to the square of the depth of color medium traversed. Such a mechanism would strongly degrade the energy of leading partons resulting in a reduced transverse momentum of leading particles in the jets that emerge after fragmentation into hadrons. The STAR experiment has shown that the topology of high pt hadron emission is consistent with jet emission, so that we may really speak about jet-suppression.

The two upper rows of fig. 9 show our  measurements of the so-called nuclear modification factors for unidentified charged hadrons from Au+Au collisions at rapidities $`\eta =0`$ and $`2.2`$ (this is the ph.d. work of Claus E. Jørgensen, see also his contribution at this conference). The nuclear modification factor is defined as:

$$R_{AA}=\frac{d^2N^{AA}/dp_td\eta }{N_{bin}d^2N^{NN}/dp_td\eta }$$
(6)
It involves a scaling of measured nucleon-nucleon transverse momentum distributions by the number of expected incoherent binary collisions, $`N_{bin}`$ (see . In the absence of any modification resulting from the ’embedding’ of elementary collisions in a nuclear collision we expect $`R_{AA}=1`$ a high $`p_t`$. At low $`p_t`$, where the particle production follows a scaling with the number of participants, the above definition of $`R_{AA}`$ leads to $`R_{AA}<1`$ for $`p_t<2GeV/c`$. In fact, it is found that $`R_{AA}>1`$ for $`p_t>2GeV/c`$ in nuclear reactions at lower energy. This effect, called the Cronin effect, is associated with initial multiple scattering of high momentum partons.

Figure 9 demonstrates that, surprisingly, $`R_{AA}<1`$ also at high $`p_t`$ for central collisions at both pseudorapidities, while $`R_{AA}1`$ for more peripheral collisions. It is remarkable that the suppression that is observed at $`p_t4GeV/c`$ is very large, amounting to a factor of 3 for central Au+Au collisions as compared to $`p+p`$ and a factor of more than 4 as compared to the more peripheral collisions. Such large suppression factors are observed at both pseudorapidities.

It has been conjectured that the observed high $`p_t`$ suppression might be the result of an entrance channel effect, for example due to a limitation in the phase space available for parton collisions related to saturation effects  in the gluon distributions inside the swiftly moving colliding nucleons (which have $`\gamma =100`$). As a test of these ideas we have very recently determined the nuclear modification factor for 100 AGeV d \+ 100 AGeV Au minimum bias collisions. The resulting $`R_{dAu}`$ is shown in fig. 10 where it is also compared to the $`R_{AuAu}`$ for central collisions previously shown in fig. 9. No high-$`p_t`$ jet suppression is observed for d+Au. In fact, the $`R_{dAu}`$ distribution measured for d+Au shows the Cronin type enhancement  observed at lower energies . At $`p_t4GeV/c`$ we find a ratio $`R_{dAu}/R_{AuAu}5`$. These observations are consistent with the smaller transverse dimensions of the overlap disk between the d and the Au nuclei and also appear to rule out strong entrance channel effects.

The very large suppression observed in central Au+Au collisions must be quantitatively understood and will require systematic dynamic modelling. At $`\eta =0`$ the particles are emitted at 90 degrees relative to the beam direction, while at $`\eta =2.2`$ the angle is only about 12 degrees. In a naive geometrical picture of an absorbing medium with cylindrical symmetry around the beam direction, the large suppression seen at forward angles indicates that the suppressing medium is extended also in the longitudinal direction. Since the observed high $`p_t`$ suppression is similar or even larger at forward rapidity as compared to midrapidity (see fig. 11) one might be tempted to infer a longitudinal extend of the dense medium which is approximately similar to its transverse dimensions ($`R5fm`$), and from this a life time longer than $`5fm/c`$. However, the problem is more complicated, due to the significant transverse and in particular longitudinal expansion that occurs as the leading parton propagates through the medium, effectively reducing the densities of color charges seen. There is little doubt that systematic studies of the high $`p_t`$-jet energy loss as a function of the thickness of the absorbing medium obtained by varying the angle of observation of high $`p_t`$ jets relative to the beam direction will be required in order to understand the properties of the dense medium. The BRAHMS experiment is uniquely suited among the RHIC experiments to carry out such a QGP tomography!

## 6 Conclusions and perspectives

The first round of RHIC experiments has clearly shown that we have moved in high energy nucleus-nucleus collisions into a qualitatively new physics domain characterized by a high degree of reaction transparency leading to the formation of a near baryon free central region. There is, nevertheless, appreciable energy loss of the colliding nuclei, so the conditions for the formation of a very high energy density zone in an interval of about $`|dN/dy|<2\pm 0.5`$ around midrapidity are present. The observation of a very large suppression of high momentum particles (jets) originating from hard scatterings in the very early stages of the collision is consistent with a significant energy loss of high momentum partons moving in a dense medium of unscreened color charges, which may be the quark gluon plasma. The next round of experiments at RHIC will afford us the opportunity to go deeper in determining the properties of the new high density state of matter that is formed, and disentangling its partonic and hadronic components . In particular it is important to determine and understand the flavor dependence of high $`p_t`$ suppression by measuring identified high $`p_t`$ mesons, and likewise to compare the suppression of baryonic and mesonic jets. A systematic scan of the jet suppression as a function of the length of dense matter traversed can be carried out in both the transverse and longitudinal directions by varying the system size and the angle of observation (rapidity). Likewise, a central future program item must be to determine, through a measurement of an excitation function, the onset of the jet suppression, which is conspicuously absent at SPS but strongly present at RHIC.

## 7 Acknowledgements

This work was supported by the Danish Natural Science Research Council, the division of Nuclear Physics of the Office of Science of the U.S. DOE, the Research Council of Norway, the Polish State Committee for Scientific Research and the Romanian Ministry of Research.


text:
# First-Principles Computations for KTN

Solid solutions of KTaO3 and KNbO3 (KTN) exhibit unusual dielectric properties at ferroelectric phase transitions, which are usually associated with a Nb instability Kleemann  that is often explained by the large Nb dynamical charge in oxides. Using non-self-consistent VASP computations Washington , we showed earlier that the Nb dynamical charge in layerd KTN strongly depends on the local environment of Nb. In the present study, we have performed self-consistent first-principles computations for fixed composition KTa7/8Nb1/8O3 in different configurations shown in Fig. 1. We computed the equilibrium ionic positions, Born charges, force constants, vibration frequencies and modes.

In agreement with our previous results, we found that the Nb dynamical charge is sensitive to the nearest neighbors (Table 1). The Nb charge in directions where there are only Nb ions is about 9.9e. In the direction where there are Ta ions, the Nb dynamical charge is about 8.3e. In contrast, the Ta charge increases in directions where there are Nb ions. The comparison of the energies (Table 2) shows that there is tendency to Nb ordering in planes if one does not take into account external stress necessary for obtaining experimental lattice constant in KTaO3. When we applied this stress, we obtained that the energies of the considered structures are close to each other. Note that this situation is very different from that known for PST, PMN etc.: in those cases, 1:1 or 1:2 order appears as a result of electrostatic interactions due to different charges of Nb and Mg or Sc and Ta Burton . In KTN, the valence of Nb and Ta is the same.

We obtained that the average Nb-O distance is larger than the average Ta-O distance. It implies that the Nb ions cause tensile distortion in KTN. Nevertheless, the lowest soft TO modes are mostly due to Nb (Table 3).

We also found that the lattice instability in the KTN 181 configuration appears in the directions of the Nb planes. That is consistent with other first-principles computations Sepliarsky . The instabilities can be explained by the long-range dipole-dipole interaction supplemented by the covalent (pseudo vibronic) effect: the displacement of one of the Nb ions along the Nb-O-Nb bond triggers appearing an additional force on nearest Nb ions in the same chain due to the change of electronic structure that is produced by the shift of the Nb ion Harrison : $`\epsilon _{cor}(R_1,R_2)=\underset{\tau ,\tau ^{}}{}\left[W(R_1)_{\tau ,\tau ^{}}W(R_2)_{\tau ^{},\tau }\right]f^\tau \left(1f^\tau ^{}\right)/\left[\epsilon ^\tau \epsilon ^\tau ^{}\right]`$ where $`W(R)=r\left(H/R\right)`$ is the linear vibronic operator, $`r`$ being the ionic displacement and $`H`$ being the electron Hamiltonian, $`\epsilon ^\tau ,f^\tau `$ are electron energy and occupation number, respectively, $`\tau `$ denotes the electronic state.

By using a tight-binding approach Nerodo ; Harrison  we have computed the correlation energy for the ionic displacements in KNbO3, in $`(z^2/R_{\mathrm{Nb}\mathrm{O}}^2)\mathrm{eV}`$ units, where $`z`$ being the ionic charge and $`R=R_{NbO}`$ is the Nb-O bond length. When the nearest Nb ions are shifted, the correlation energy is –10.3 and –0.9 for the displacements along and perpendicular to the Nb-O-Nb line, respectively. If Nb and nearest O are shifted then the correlation energy is 25.9 and 8.6 for the displacemnts along and perpendicular to the Nb-O bond, respectively. Hence, the vibronic effect facilitates the displacement of the Nb ions in the chains in the same direction while the Nb and O ions are easily shifted in opposite directions. This result is in qualitative agreement with straightforward computations of the interatomic force constants (see Tinte  and references therein), which showed that the M-O force constant along the …-M-O-M-O-…chains proves to be practically negligible while the M-M force constant appears to be surprisingly strong.

Hence, the local polarization in KTN appears in the Nb-reached regions. The polarization direction follows the field and the free energy can be written as: $`F=\alpha (PP_0)^2+A(PP_0)^3+\mathrm{}(E+E_0)P`$ where $`E_0`$ is random and $`P`$ is polarization *magnitude* ($`P>P_0>0`$ where $`P_0`$ is the lowest possible, remnant, magnitude of polarization). An unusual field dependence of susceptibility follows from this: $`1/\chi (E)=2\alpha +3AE/\alpha +\mathrm{},EE_0`$. At smaller fields, susceptibility is constant or decreases because of averaging over $`E_0`$ hydro .

The strongly angle dependent correlation of the displacements can explain the diffuse X-ray scattering in KNbO3 Comes  without the assumption that the Nb ions have deep off-center wells. At high temperatures (where the soft mode is stabilized), acoustic excitations of the (Nb or O) chains can be intensive due to nearly flat transverse acoustic phonon energy surface in (001) planes Dorner . The peaks in the XANES data Ravel  can be understood if the large chain-like fluctuations are considered together with the strong dependence of the X-Ray absorption matrix element on the ionic displacements Vedrinskii .

The dynamical nature of the cooperative displacements Takahasi  is supported by the fact that, in the high temperature cubic phase, there is little dielectric dispersion Turik . In the ferroelectric phases, the situation is more complicated. For instance, in the tetragonal phase, the average polarization is directed along the  direction. However, the direction of a local Nb displacement is  due to anharmonicity. The correlation of the displacements can result in the appearance of finite size chains directed in the  directions, while the metal ions are displaced in  directions (common within a finite-size chain but random in different chains Takahasi  ). One can introduce population numbers for such chains corresponding to different directions of the displacements. Note that these (quasistatic) chains can exist only because of i) anharmonicity of the potential relief, where the direction of the (local) displacement does not coincide with the direction of the polarization, ii) strong correlation of the displacements, iii) existence of a uniform macroscopic field. In the cubic phase, the last condition fails and the dispersion connected with the chains disappears. In PbTiO3, the low-temperature phase is tetragonal: the ions are shifted in the  direction and the first condition fails; the dispersion is also absent in accordance with experiment Turik . In the ferroelectric phases of BaTiO3 and KNbO3, all three conditions are satisfied and the dispersion exists Turik .


text:
# Geo-neutrinos, Mantle Circulation and Silicate Earth

## 1 Introduction

KamLAND, the Low Energy Anti-Neutrino Detector operating at the Kamioka mine in Japan, has demonstrated the possibility  of detecting geo-neutrinos, the antineutrinos from decay chains of radiogenic nuclides inside Earth. A new window on Earth’s interior is thus being opened, see, *e.g.*, Refs. .

Recently, a reference model of geo-neutrino fluxes has been presented in Ref. . The model, hereafter indicated as REF, is based on a detailed description of Earth’s crust and mantle, and takes into account available information on the abundances of Uranium, Thorium and Potassium — the most important heat and neutrino sources — inside Earth’s layers.

This model has to be intended as a starting point, providing first estimates of expected events at several locations on the globe, see Fig. 1 111An interactive version of the geo-neutrino event map shown in Fig. 1 will be available at the site http://www.fe.infn.it/$`\stackrel{~}{}`$mantovani/geoneutrini/index.htm .. In preparation to the experimental results which will be available in the future, from KamLAND as well as from other detectors which are in preparation, it is useful to consider geo-neutrino production in greater detail, putting the basis for a more refined model. In this respect, let us discuss critically the key ingredients of REF.

Concerning the crust, the $`2^{}\times 2^{}`$ model of Ref.  was adopted for building REF; world-average abundances $`a`$ of radiogenic elements have been estimated separately for oceans, the continental crust (subdivided into upper, middle and lower sub-layers), sediments, and oceanic crust. Although this treatment looks rather detailed on the globe scale, the typical linear dimension of each tile is of order 200 km, so that any information on a smaller scale is essentially lost. Since one expects that about one third of the signal is contributed from sources within 100 Km from the detector, a better description of the surrounding region is needed for a refined estimate.

Sandwiched between Earth’s crust and metallic core, the mantle is a 2900 km layer of pressurized rock at high temperature. As reviewed by Hofmann in Refs. , mantle models can be divided into two broad classes, essentially corresponding to the presently contradictory geochemical and geophysical evidence of Earth’s interior.

Geochemists have long insisted on a two-layer model, in which the mantle consists of a relatively primitive layer below a depth of about 670 km and an upper layer that is highly depleted of heat producing elements (panel “a” in Fig. 2). The two layers are viewed as separate sources of the Mid-Ocean-Ridge Basalts (MORB), which come from mantle regions that have been already depleted in incompatible elements by extraction of the continental crust, and of Ocean Island Basalts (OIB), which form by melting of deeper, less depleted or even enriched mantle sectors. Also, a more primitive deep layer is needed from global constraints, otherwise the amount of radiogenic elements inside Earth is much too small with respect to that estimated within the Bulk Silicate Earth (BSE) paradigm (see below).

On the other hand, over the past several years seismic tomography has provided increasingly detailed images of apparently cold slab descending into the deep mantle, below the 670-km boundary. If cold slabs descend into the deep mantle, there must be a corresponding upward flow of deep-mantle material to shallow levels (panel “b” in Fig. 2). If this circulation reaches the bottom of the mantle (whole mantle convection), it would destroy any compositional layering below the crust in a few hundred million years (at a typical speed of 3 cm yr-1 it takes about 108 yr to move down to 2900 km).

In REF, the Preliminary Reference Earth Model  was used for the mantle density profile, dividing Earth’s interior into several spherically symmetrical shells corresponding to seismic discontinuities. Concerning its composition, a two-layer stratified model was used: for present day upper mantle, considered as the source of MORB, mass abundances of 6.5 and 17.3 ppb for Uranium and Thorium, respectively, and 78 ppm for Potassium were assumed down to a depth $`h_0=670`$ km. These abundances were obtained by averaging the results of Refs.  and . Abundances below were inferred by requiring that the BSE constraint is globally satisfied, obtaining 13.2 and 52 ppb for U and Th, 160 ppm for K. In view of the present debate about mantle circulation and composition, it is desirable to have a more general treatment, which encompasses both geochemical and geophysical supported models.

The Bulk Silicate Earth framework provides a compositional description of the primitive mantle, the outer part of Earth after core separation and before crust differentiation. In this frame, the amount of heat/neutrino generating material inside Earth is determined through the following steps:

(a) From the compositional study of selected samples emerging from the mantle, after correcting for the effects of partial melting, one establishes the absolute primitive abundances in major elements with refractory and lithophile character, *i.e.*, elements with high condensation temperature (so that they do not escape in the processes leading to Earth formation) and which do not enter the metallic core. In this way primitive absolute abundances of elements such as Al, Ca and Ti are determined, a factor about 2.75 times CI chondritic abundances.

(b) It is believed, and supported by studies of mantle samples, that refractory lithophile elements inside Earth are in the same proportion as in chondritic meteorites. In this way, primitive abundances of Th and U can be derived by rescaling the chondritic values 222Note that Uranium absolute abundances in chondritic meteorites are variable within a factor two, however the ratio to Aluminum is stable within 10%..

(c) Potassium, being a moderately volatile elements, could have escaped in the planetesimal formation phase. Its absolute abundance is best derived from the practically constant mass ratio with respect to Uranium observed in crustal and mantle derived rocks.

The BSE estimates for the total amounts of Uranium, Thorium and Potassium from different authors  are quite concordant within 10%, the central values being $`m_{\mathrm{BSE}}=0.8\times 10^{17}`$ kg for Uranium, $`3.1\times 10^{17}`$ kg for Thorium, and $`0.9\times 10^{21}`$ kg for Potassium. These values can be taken — within their uncertainties — as representative of the composition of the present crust plus mantle system 333 Geochemists generally agree on the absence of these lithophile elements in Earth’s core, see however Refs.  for a different point of view..

Different models can provide different distributions between crust and mantle, however for each element the sum of the masses is fixed by the BSE constraint. This clearly provides constraints on the geo-neutrino flux 444We are always speaking of the flux integrated over any direction, see  for a precise definition of the geo-neutrino flux. which are grounded on sound geochemical arguments or, alternatively, geo-neutrino detection can provide a test of an important geochemical paradigm.

We remark that BSE predictions for the radiogenic contributions to Earth’s flow are $`H(\mathrm{U})=7.6`$ TW , $`H(\mathrm{Th})=8.5`$ TW and $`H(\mathrm{K})=1.8`$ TW, their sum being about one half of the observed heat flow from Earth, $`H_E40`$ TW. Geo-neutrino detection can thus provide a direct insight on the main source of Earth’s energetics, see .

The reference model was built so as to satisfy the BSE constraint. However, when discussing its uncertainties a wide range of models was considered, from a minimal model providing a radiogenic contribution of just 9 TW up to a fully radiogenic model. These extreme models, although not excluded from direct observational data, are well outside the BSE geochemical constraint. If the BSE global constraint is imposed and only the crust contribution is left free within the observational uncertainties, then the geo-neutrino prediction becomes much more precise.

Briefly, in this paper we shall address three questions.

(i) How sensitive are the predicted geo-neutrino fluxes to uncertainties about the mechanism of mantle circulation?

(ii) Is it possible to test the Bulk Silicate Earth model with neutrinos?

(iii) Which accuracy is needed concerning the local/regional predictions and the geo-neutrino detection?

We shall restrict the discussion to geo-neutrinos from Uranium progeny, which are more easily detectable due to the higher energy. Extension to the other chains is immediate.

## 2 Geochemistry, geophysics and geo-neutrinos

As mentioned in the introduction, the composition and circulation inside Earth’s mantle is the subject of a strong and so far unresolved debate between geochemists and geophysicists. Geochemistry supports the existence of two compositionally distinct reservoirs in the mantle, the borders between them being usually placed at depth near $`h_0=670`$ km, whereas geophysics presents evidence of matter circulation extending well beyond this depth. If this circulation involves the whole mantle, it would have destroyed any pre-existing layering, in conflict with geochemical evidence. In this section we look at the implications of this controversy on the predicted geo-neutrino fluxes.

One can build a wide class of models, including the extreme geochemical and geophysical models, in terms of just one free parameter, the depth $`h`$ marking the borders between the two hypothetical reservoirs:

(i) We assume that Uranium abundance in the uppermost part of the mantle is nearly chondritic 555We are aware that current estimates of U in depleted upper mantle after crust extraction are in the range of 2 to 3.9 ppb . Given the uncertainty on these values we prefer to perform calculations using the well-constrained CI value. As shown below, the assumption of lower U abundance for the uppermost depleted mantle has limited effects on geo-neutrino flux predictions.  ($`a_u=6.5`$ ppb) down to an unspecified depth $`h`$.

(ii) Below $`h`$ we determine abundances ($`a_l`$) by requiring mass balance for the whole Earth.

This means that Uranium mass below the critical depth, $`m_{>h}`$, is obtained by subtracting from the total BSE estimated mass ($`m_{\mathrm{BSE}}`$) the quantity observationally determined in the crust ($`m_c`$) and that contained in the mantle above $`h`$ ($`m_{<h}`$):

$$m_{>h}=m_{\mathrm{BSE}}m_cm_{<h}$$
(1)
The abundance in the lower part is then calculated as the ratio of $`m_{>h}`$ to Earth’s mass below $`h`$ ($`M_{>h}`$):

$$a_l=m_{>h}/M_{>h}$$
(2)
This class of models, described in Figs. 4 and 4, includes a fully mixed mantle, which is obtained for $`h=25`$ km (*i.e.*, just below the crust) so that the strongly impoverished mantle has a negligible thickness. The traditional geochemical model corresponds to $`h=h_0`$. As $`h`$ increases, the depleted region extends deeper inside the Earth and — due to mass balance — the innermost part of the mantle becomes richer and closer in composition to the primitive mantle.

Let us discuss in detail a few cases, remembering that the BSE estimate for Uranium in the whole Earth is $`m_{\mathrm{BSE}}=0.8\times 10^{17}`$ kg and that the best estimate for the amount in the crust is $`m_c=0.35\times 10^{17}`$ kg  so that Uranium in the mantle is $`m_m=0.45\times 10^{17}`$ kg.

(a) In the fully mixed model, this quantity has to be distributed over the mantle mass $`M_m=4.0\times 10^{24}`$ kg, which yields a uniform mantle abundance $`a=11.25`$ ppb. We shall refer to this model as MIX.

(b) If we keep the estimated abundance in the uppermost part ($`a_u=6.5`$ ppb) down to $`h_0`$ one has the REF model .

(c) Among all possible models, the case $`h=1630`$ km is particularly interesting. Below this depth the resulting Uranium abundance is 20 ppb, corresponding to the BSE estimate. The innermost part of the mantle is thus primitive in its trace element composition and the crust enrichment is obtained at expenses of the mantle content above $`h`$. We shall refer to this model as PRIM.

Concerning geo-neutrino fluxes from the mantle, all these models have the same amount of heat/neutrino sources and only the geometrical distribution is varied. The largest flux corresponds to the model with more sources close to the surface, *i.e.* to the MIX model. On the other hand, the minimal prediction is obtained when the sources are most concentrated at larger depth, which corresponds to the PRIM case.

From Table 1, the difference between the extreme cases is 8%, model REF being in between. The abundance in the upper reservoir $`a_u`$ can also be treated as a free parameter. If we use an extremely low value $`a_u=2`$ ppb  down to about 1200 km and primitive abundance below, we obtain the minimal prediction $`0.86\times 10^6`$ cm-2 s-1.

We conclude this section with the following remarks:

* uncertainties on the geometrical distribution of trace elements in the mantle can change the REF prediction for the mantle by at most $`\pm 8\%`$.
* A geo-neutrino detector at a site where the contribution from the mantle is dominant (*i.e.*, far from the continental crust) can be sensitive to the mantle compositional geometry only if a percent accuracy is reached.
* Since at Kamioka mine or at Gran Sasso the mantle contribution to the total flux is about one quarter of the total , uncertainties on the mantle geometry imply an estimated error of about 2% on the total flux predicted with REF.

## 3 The Bulk Silicate Earth constraint

So far we have been considering the effect of different geometrical distributions of trace elements in the mantle, for fixed amount of these elements within it. Actually the BSE model can be exploited to obtain tight constraints on the total flux produced together from the crust and the mantle. In fact, with BSE fixing the total amount of trace elements inside Earth, geometrical arguments and observational constraints on the crust composition can be used in order to find extreme values of the produced fluxes. As an extension of the previous section, the maximal (minimal) flux is obtained by placing the sources as close (far) as possible to Earth’s surface, where the detector is located.

Concerning Uranium, the range of BSE concentrations reported in the literature is between 18 and 23 ppb, corresponding to a total Uranium mass between $`m(\mathrm{min})=0.72`$ and $`m(\mathrm{max})=0.92`$ in units of $`10^{17}`$ kg.

In the same units, we estimate that Uranium mass in the crust is between $`m_c(\mathrm{min})=0.30`$ and $`m_c(\mathrm{min})=0.41`$, by taking the lowest (highest) concentration reported in the literature for each layer, see Table II of . The main source of uncertainty is from the abundance in the lower crust, estimated at 0.20 ppm in Ref.  and at 1.1 ppm in Ref.  Estimates for the abundance in the upper crust are more concordant, ranging from 2.2 ppm  to 2.8 ppm . We remark that, within this appraoch, the resulting *average* crustal U abundance $`a_{cc}`$ is in the range 1.3-1.8 ppm, which encompasses all estimates reported in the literature  but for that of Ref. , $`a_{cc}=0.91`$ ppm, see Table 2 666Note that crust mass and flux ranges found in this paper are tighter than those quoted in Ref. , which used the value $`a_{cc}=0.91`$ ppm from Ref.  as lower limit..

The highest flux is obtained by assuming the maximal mass in the crust and the maximal mass in the mantle, $`m_{\mathrm{tot}}(\mathrm{max})m_c(\mathrm{max})=0.51`$, with a uniform distribution inside the mantle, corresponding to $`a=12.8`$ ppb.

On the other hand, the lowest flux corresponds to the minimal mass in the crust and the minimal mass in the mantle, $`m_{\mathrm{tot}}(\mathrm{min})m_c(\mathrm{min})=0.42`$, with a distribution in the mantle similar to that of PRIM, *i.e.* a strongly depleted mantle with $`a_u=2`$ ppb down to about 1300 km and a primordial composition beneath.

The predicted fluxes are shown in Table 3 for a few locations of particular interest: the Kamioka mine ($`33^{}`$N $`85^{}`$E) where KamLAND is operational, the Gran Sasso laboratory ($`42^{}`$N $`14^{}`$E) where BOREXINO  is in preparation, the top of Himalaya ($`36^{}`$N $`137^{}`$E), which receives the maximal contribution from the crust, and Hawaii ($`20^{}`$N $`156^{}`$E), a location where the mantle contribution is dominant.

At any site the difference between the maximal and the minimal flux predictions are of about 30%, the extreme values being within $`\pm 15\%`$ from the reference model prediction.

All this shows the power of the BSE constraint. If the total amount of Uranium inside Earth is fixed at $`m_{\mathrm{BSE}}=(0.8\pm 0.1)\times 10^{17}`$ kg, then the produced geo-neutrino flux at, *e.g.*, Kamioka is

$$\mathrm{\Phi }=(3.7\pm 0.6)\times 10^6\mathrm{cm}^2\mathrm{s}^1(\mathrm{full}\mathrm{range})$$
(3)
after taking into account the full range of global observational uncertainties on Uranium abundances in the crust and uncertainties concerning circulation in the mantle.

For a comparison, we remind that in Ref. , where the Uranium mass was allowed to vary in the range (0.4–1.7)$`\times 10^{17}`$ kg, the whole range uncertainty on the produced flux at Kamioka was $`\pm 1.9\times 10^6`$ cm-2s-1.

We remark that the error quoted in Eq. (3) corresponds to a full range of the predicted values. If, following a commonly used rule of thumb, we consider such a range as a $`\pm 3\sigma `$ (99.5%) confidence level, we deduce a conventional $`1\sigma `$ estimate:

$$\mathrm{\Phi }=(3.7\pm 0.2)\times 10^6\mathrm{cm}^2\mathrm{s}^1(1\sigma ).$$
(4)
## 4 The effects of geochemical fluctuations and neutrino oscillations

The main result of the previous section is that — neglecting regional fluctuations — global mass balance provides a precise determination of the produced geo-neutrino fluxes. We shall compare this precision with uncertainties resulting from fluctuations of the regional geochemical composition, from the available information on neutrino mixing parameters and from the detector finite size.

(a) Actually the Uranium concentration in the region where the detector is located can be different from the world average and local fluctuations of this highly mobile element can be envisaged. These variations, although negligible for mass balance, can affect significantly the flux. In other words, geochemical arguments fix the contribution of distant sources, and a more detailed geological and geochemical investigation of the region around the detector is needed, the error quoted in Eq. (4) providing a benchmark for the accuracy of the local evaluation.

As an example, it has been estimated that about one third (one half) of the geo-neutrino signal is generated within a distance of 100 (500) km from Kamioka, essentially in the Japanese continental shelf. In REF the world averaged upper crust Uranium concentration $`a_{uc}=2.5`$ ppm was adopted for Japan. In a recent study of the chemical composition of Japan upper crust  more than hundred samples, corresponding to 37 geologic groups, have been analyzed. The composition is weighted with the frequency in the geological map and the resulting average abundance is $`a_{\mathrm{Jap}}=2.32`$ ppm , which implies a 7.2% reduction of the flux from Japanese upper crust with respect to the estimate in REF. Larger variations occur when rocks are divided according to age or type, see Table 4, and even larger differences are found within each group. All this calls for a detailed geochemical and geophysical study, with the goal of reducing the effect of regional fluctuations to the level of the uncertainty from global geochemical constraints.

(b) Let us remark that the signal is originated from neutrinos which maintain the electron flavour in their trip from source to detector, the effective flux being $`\mathrm{\Phi }_{\mathrm{eff}}=\mathrm{\Phi }P_{ee}`$, where $`\mathrm{\Phi }`$ is the produced flux and $`P_{ee}=11/2\mathrm{sin}^2(2\theta )`$ is the (distance averaged) survival probability. From the analysis of all available solar and reactor neutrino experiments, one gets $`\mathrm{tan}^2\theta =0.41\pm 0.04`$ at $`1\sigma `$ , which implies:

$$P_{ee}=0.59\pm 0.02(1\sigma ).$$
(5)
The survival probability is presently know with a 3% accuracy, so that uncertainties on the neutrino fate do not mask its sources.

(c) If Uranium geo-neutrinos are detected by means of inverse beta reaction on free hydrogen nuclei ($`\overline{\nu }_e+pe^++n`$) the event number is

$$N=13.2ϵ\frac{\mathrm{\Phi }_{\mathrm{eff}}}{10^6\mathrm{cm}^2\mathrm{s}^1}\frac{N_p}{10^{32}}\frac{t}{\mathrm{yr}},$$
(6)
where $`ϵ`$ is the detection efficiency, $`N_p`$ is the number of free protons in the target and $`t`$ is the measurement time. For a produced flux $`\mathrm{\Phi }=4\times 10^6`$ cm-2 s-1 and $`ϵ=0.8`$, one expects 25 events for an exposure of $`10^{32}`$ protons yr. Statistical fluctuations will be of order $`\sqrt{N}`$ if background can be neglected (this clearly does not hold for Kamioka, due to the many nearby nuclear power plants). In order to reach a 5% accuracy — comparable to that of the global geochemical estimate — one needs an exposure of $`16\times 10^{32}`$ protons yr, which corresponds to a five-kton detector operating over four years.

## 5 Concluding remarks

We summarize here the main points of this paper:

(1) uncertainties on the geometrical distribution of trace elements in the mantle (for a fixed mass within it) can change the prediction of the reference model  for the geo-neutrino flux from mantle by at most $`\pm 8\%`$ (full range), the extreme values corresponding to a fully-mixed and to a two-layer model, with primordial abundance below about 1300 km.

(2) By using global mass balance for the Bulk Silicate Earth, the predicted flux contribution originating from distant sources in the crust and in the mantle is fixed within $`\pm 5\%`$ ($`1\sigma `$) with respect to the reference model.

(3) A detailed geological and geochemical investigation of the region within few hundreds km from the detector has to be performed, for reducing the flux uncertainty from fluctuations of the local abundances to the level of the global geochemical error.

(4) A five-kton detector operating over four years at a site relatively far from nuclear power plants can measure the geo-neutrino signal with 5% accuracy.

(5) This will provide a crucial test of the Bulk Silicate Earth and a direct estimate of the radiogenic contribution to terrestrial heat.

###### Acknowledgments.

We express our gratitude for useful discussions Dr. C. Bonadiman, L. Carmignani, M. Coltorti, S. Enomoto, K. Inoue, E. Lisi, T. Mitsui, B. Ricci, N. Sleep, A. Suzuki, and F. Villante. G. F. is grateful to the organizing committee of AHEP 2003 for a very enjoiable conference. This work was partially supported by MIUR (Ministero dell’Istruzione, dell’Università e della Ricerca) under MIUR-PRIN-2003 project “Theoretical Physics of the Nucleus and the Many-Body Systems” and MIUR-PRIN-2002 project “Astroparticle Physics”.