format: text

text:
# Defining and identifying communities in networks

## I Introduction

Evidence has rapidly grown in the past few years that very diverse systems in many different fields can be described as complex networks, i. e. assemblies of nodes and edges with nontrivial topological properties Barabasi02 ; Newman03 . Examples range from technological systems (the Internet and the web AJB99 ; Broder00 ) to biological (epidemiology MN00 ; PV01 , metabolic networks Jeong00 ; WF01 ; barabasi\_mod , ecological webs Montoya02 ; DWM02 ; CGA02 ; gcalda ) and social systems Redner ; Newman01bis  (scientific collaborations, structure of large organizations).

In this paper we deal with a topological property of networks, the community structure, that has attracted a great deal of interest very recently. The concept of community is very common and it is linked to the classification of objects in categories for the sake of memorization or retrieval of information. From this point of view the notion of community is very general and, depending on the context, can be synonymous of module, class, group, cluster, etc. Among the many contexts where this notion is relevant it is worth mentioning the problem of modularity in metabolic or cellular networks barabasi\_mod ; rives  or the problem of the identification of communities in the web webcomm . This last issue is relevant for the implementation of search engines of new generation, content filtering, automatic classification or the automatic realization of ontologies.

Given the relevance of the problem it is crucial to construct efficient procedures and algorithms for the identification of the community structure in a generic network. This, however, is a highly nontrivial task.

Qualitatively, a community is defined as a subset of nodes within the graph such that connections between the nodes are denser than connections with the rest of the network. The detection of the community structure in a network is generally intended as a procedure for mapping the network into a tree (Fig. 1). In this tree (called dendrogram in social sciences) the leaves are the nodes while the branches join nodes or (at higher level) groups of nodes, thus identifying a hierarchical structure of communities nested within each other.

Several algorithms to perform this mapping are known in the literature. The traditional method is the so-called hierarchical clustering Wassermanbook . For every pair $`i,j`$ of nodes in the network one calculates a weight $`W_{i,j}`$, which measures how closely connected the vertices are. Starting from the set of all nodes and no edges, links are iteratively added between pairs of nodes in order of decreasing weight. In this way nodes are grouped into larger and larger communities and the tree is built up to the root, which represents the whole network. Algorithms of this kind are called agglomerative.

For the other class of algorithms, called divisive, the order of construction of the tree is reversed: one starts with the whole graph and iteratively cuts the edges, thus dividing the network progressively into smaller and smaller disconnected sub-networks identified as the communities. The crucial point in a divisive algorithm is the selection of the edges to be cut, which have to be those connecting communities and not those within them. Very recently, Girvan and Newman (GN) have introduced a divisive algorithm where the selection of the edges to be cut is based on the value of their “edge betweenness” Girvan02 , a generalization of the centrality betweenness introduced by Anthonisse Anthonisse71  and Freeman Freeman77 . Consider the shortest paths between all pairs of nodes in a network. The betweenness of an edge is the number of these paths running through it. It is clear that when a graph is made of tightly bound clusters, loosely interconnected, all shortest paths between nodes in different clusters have to go through the few inter-clusters connections, which therefore have a large betweenness value. The single step of the GN detection algorithm consists in the computation of the edge betweenness for all edges in the graph and in the removal of those with the highest score. The iteration of this procedure leads to the splitting of the network into disconnected subgraphs that in their turn undergo the same procedure, until the whole graph is divided in a set of isolated nodes. In this way the dendrogram is built, from the root to the leaves.

The GN algorithm represents a major step forward for the detection of communities in networks, since it avoids many of the shortcomings of traditional methods Girvan02 ; Newman\_cond-mat . This explains why it has been quickly adopted in the past year as a sort of standard for the analysis of community structure in networks holme ; Wilkinson02 ; gleiser ; guimera ; tyler .

This paper follows a different track by proposing an alternative strategy for the identification of the community structure. This complementary approach follows from the need of addressing the two following issues.

1. In general algorithms define communities operationally as what the they finds. A dendrogram, i. e. a community structure, is always produced by the algorithms down to the level of single nodes, independently from the type of graph analyzed. This is due to the lack of explicit prescriptions to discriminate between networks that are actually endowed with a community structure and those that are not. As a consequence, in practical applications one needs additional, non topological, information on the nature of the network to understand which of the branches of the tree have a real significance. Without such information it is not clear at all whether the identification of a community is reliable or not. There have been two noticeable proposals to solve this problem. In particular it is worth mentioning the approach proposed by Wilkinson and Huberman Wilkinson02 , which is limited to the lowest level of the community structure and specific to algorithms based on betweenness. More recently Newman\_cond-mat , Newman and Girvan have introduced an a posteriori measure of the strength of the community structure, which they called modularity. More precisely, the modularity estimates the fraction of inward links in a community minus the expectation value of the same quantity in a network with the same community divisions but random connections between the nodes. This quantity definitely gives an indication of the strength of the community structure, even though the lack of the implementation of a quantitative definition of community does not allow to discriminate in an objective way meaningful communities.
2. The “edge betweenness algorithm” is computationally costly, as already remarked by Girvan and Newman Girvan02 ; Newman\_cond-mat . Evaluating the score for all edges requires a time of the order of $`MN`$, where $`M`$ is the number of edges and $`N`$ the number of nodes. The iteration of the procedure for all $`M`$ edges leads in the worst case to a total scaling of the computational time as $`M^2N`$, which makes the analysis practically unfeasible already for moderately large networks (of the order of $`N=10000`$ Newman\_cond-mat ).

In this paper we propose solutions to both these problems. First we introduce a general criterion for deciding which of the subgraphs singled out by the detection algorithms are actual communities. We discuss in detail the case of two quantitative definitions of community. In this way we transform the GN algorithm in a self-contained tool. Secondly we present an alternative algorithm, based on the computation of local quantities, which gives, in controlled cases, results of accuracy comparable to the GN method, while largely outperforming it from the point of view of the computational speed.111It is worth mentioning that after the completion of this work Newman has proposed a new agglomerative algorithm to address the issue of the computational efficiency newman\_fast .

The outline of the paper is as follows. In section II we introduce the definitions of community. In section III we show how these definitions can be implemented in a generic divisive algorithm in order to make it self-contained and we present tests on some computer-generated and real networks. In section IV we present a new and fast algorithm for the detection of communities and we compare its performance with the GN algorithm. Section V is devoted to the application of the new algorithm to a network of scientific collaborations which, because of its size, is hard to analyze with the Girvan-Newman method. We finally draw some conclusions and discuss the perspectives of this work.

## II Quantitative definitions of community

The idea to solve the first of the problems discussed above is very simple: the algorithm that builds the tree just selects subgraphs that are candidate to be considered communities. One has then to check whether they are actually such by using a precise definition. If the subgraph does not meet the criterion, the subgraph isolated from the network is not a community and the corresponding branch in the dendrogram should not be drawn.

As mentioned above, a community is generally thought as a part of a network where internal connections are denser than external ones. To sharpen the use of detection algorithms a more precise definition is needed. Many possible definitions of communities exist in the literature Wassermanbook . Here we consider explicitly the implementation in the algorithms of two plausible definitions of community which translate into formulas the sentence above.

The basic quantity to consider is $`k_i`$, the degree of a generic node $`i`$, which in terms of the adjacency matrix $`A_{i,j}`$ of the network $`G`$ is $`k_i=_jA_{i,j}`$ 222The adjacency matrix fully specifies the topology of the network. In the simplest case of an unweighted, undirected network, it is equal to 1 if $`i`$ and $`j`$ are directly connected, zero otherwise.. If we consider a subgraph $`VG`$, to which node $`i`$ belongs, we can split the total degree in two contributions: $`k_i(V)=k_i^{in}(V)+k_i^{out}(V)`$. $`k_i^{in}(V)=_{jV}A_{i,j}`$ is the number of edges connecting node $`i`$ to other nodes belonging to $`V`$. $`k_i^{out}(V)=_{jV}A_{i,j}`$ is clearly the number of connections toward nodes in the rest of the network.

Definition of community in a strong sense

The subgraph $`V`$ is a community in a strong sense if

$$k_i^{in}(V)>k_i^{out}(V),iV.$$
(1)
In a strong community each node has more connections within the community than with the rest of the graph. This definition coincides with the one proposed in webcomm  in the framework of the identification of web communities.

Definition of community in a weak sense

The subgraph $`V`$ is a community in a weak sense if

$$\underset{iV}{}k_i^{in}(V)>\underset{iV}{}k_i^{out}(V).$$
(2)
In a weak community the sum of all degrees within $`V`$ is larger than the sum of all degrees toward the rest of the network.

Clearly a community in a strong sense is also a community in a weak sense, while the converse is not true.

It is worth mentioning that our definitions of community, though very natural, do not represent the only possible choice. Several other possible definitions, possibly more appropriate in some cases, exist and are described in Wassermanbook . Among them for instance the definition of the so-called $`LSset`$ goes in the direction of our strong definition even though extremely more stringent. An $`LSset`$ is a set of nodes such that each of its proper subsets has more ties to its complement within the set than outside. On the other hand the definition of $`k`$-core is roughly, although not exactly, equivalent to our weak definition. A $`k`$-core is defined as a subgraph in which each node is adjacent to at least a minimum number, $`k`$, of the other nodes in the subgraph.

## III Self-contained algorithms

From the definitions given above, it is apparent that, if a network is randomly split in two parts, one very large and the other with only few nodes, the very large part almost always fulfills the definition of community. In order to deal with this problem, let us consider the Erdös-Renyi random graph Bollobasbook . If we cut at random the graph in two parts containing $`\alpha N`$ and $`(1\alpha )N`$ nodes, respectively, it is easy to evaluate analytically the probability $`P(\alpha )`$ that the subgraph containing $`\alpha N`$ nodes fulfills the weak or the strong definition (details will be given in a forthcoming publication). It turns out that, as soon as $`N`$ is sufficiently large, the probability is very close to a step function around $`\alpha =0.5`$. Hence it is extremely likely that, in a random graph randomly cut in two parts, the largest one is a community according to the previous definitions. However it is extremely unlikely that both subgraphs fulfill simultaneously the definitions: therefore if we accept divisions only if both groups fulfill the definition of community, we correctly find that a random graph has no community structure. We extend this criterion to the general case: if less than two subgraphs obtained from the cut satisfy the definitions, then the splitting is considered to be an artifact and disregarded.

We can now summarize the improved self-contained version of the GN algorithm.

1. Choose a definition of community;
2. Compute the edge betweenness for all edges and remove those with the highest score;
3. If the removal does not split the (sub-)graph go to point 2;
4. If the removal splits the (sub-)graph, test if at least two of the resulting subgraphs fulfill the definition. If they do, draw the corresponding part of the dendrogram;
5. Iterate the procedure (going back to point 2) for all the sub-graphs until no edges are left in the network.

It is important to remark that the quantities appearing in Eqs. (1) and (2) must always be evaluated with respect to the full adjacency matrix. The application of this procedure to a network produces a tree, where every branch splitting represents a meaningful (with respect to the definition) separation in communities.

It is now possible to blindly test the effectiveness of the GN algorithm. We have considered the “artificial” graph already discussed by Girvan and Newman. It is a simple network with $`N`$ nodes divided into four groups: connections between pairs within a group are present with probability $`p_{in}`$, while pairs of nodes in different groups are connected with probability $`p_{out}`$. As the probability $`p_{out}`$ grows from zero, the community structure in the network becomes less well defined.

For every realization of the artificial graph the application of the detection algorithm generates a tree. We consider the algorithm to be successful if the four communities are detected, each node is classified in the right community and the communities are not further subdivided.

Fig. 2 presents the comparison of the fraction of successes for the modified GN algorithm with the expected value computed analytically. We see that the GN algorithm captures very well the existence of communities in a strong sense, while it performs less well for the weak definition. However, one should not be misguided by the quantity presented in Fig. 2. By looking at a softer measure of success, the fraction $`f`$ of nodes not correctly classified, one realizes that, when the algorithm with weak definition seems to fail, it correctly identifies the four communities and it misclassifies only a few nodes up to much higher values of $`p_{out}`$. The deviations from the theoretical behaviour observed for small values of $`p_{out}`$ are due to the possibility that one or more of the four communities are further split in smaller sub-communities. This event, not taken into account in the analytical calculation, becomes very unlikely as the size of the system increases.

## IV A new fast algorithm

The Girvan-Newman algorithm is computationally expensive because it requires the repeated evaluation, for each edge in the system, of a global quantity, the betweenness, whose value depends on the properties of the whole system. Despite smart methods to compute the edge betweenness simultaneously for all edges Newman01 ; Brandes01 , the evaluation of such quantity is the time consuming part of the procedure. As a consequence the time to analyze completely a network turns out to grow fast with its size, making the analysis unfeasible for networks of size larger than around $`10000`$ nodes Newman\_cond-mat .

In order to overcome this problem we introduce a new kind of divisive algorithm which requires the consideration of local quantities only, and is therefore much faster than the GN algorithm. The fundamental ingredient of a divisive algorithm is a quantity which can single out edges connecting nodes belonging to different communities. We consider the edge clustering coefficient, defined, in analogy with the usual node clustering coefficient, as the number of triangles to which a given edge belongs, divided by the number of triangles that might potentially include it, given the degrees of the adjacent nodes. More formally, for the edge connecting node $`i`$ to node $`j`$, the edge clustering coefficient is

$$C_{i,j}^{(3)}=\frac{z_{i,j}^{(3)}}{\mathrm{min}[(k_i1),(k_j1)]},$$
(3)
where $`z_{i,j}^{(3)}`$ is the number of triangles built on that edge and $`\mathrm{min}[(k_i1),(k_j1)]`$ is the maximal possible number of them.

The idea behind the use of this quantity in a divisive algorithm is that edges connecting nodes in different communities are included in few or no triangles, and tend to have small values of $`C_{i,j}^{(3)}`$. On the other hand many triangles exist within clusters. Hence the coefficient $`C_{i,j}^{(3)}`$ is a measure of how inter-communitarian a link is. A problem arises when the number of triangles is zero, because $`C_{i,j}^{(3)}=0`$, irrespective of $`k_i`$ and $`k_j`$, (or even $`C_{i,j}^{(3)}`$ is indeterminate, when $`\mathrm{min}[(k_i1),(k_j1)]=0`$). To remove this degeneracy we consider a slightly modified quantity by using, at the numerator, the number of triangles plus one:

$$\stackrel{~}{C}_{i,j}^{(3)}=\frac{z_{i,j}^{(3)}+1}{\mathrm{min}[(k_i1),(k_j1)]}.$$
(4)
By considering higher order cycles we can define, in much the same way, coefficients of order $`g`$ as:

$$\stackrel{~}{C}_{i,j}^{}{}_{}{}^{(g)}=\frac{z_{i,j}^{(g)}+1}{s_{i,j}^{(g)}},$$
(5)
where $`z_{i,j}^{(g)}`$ is the number of cyclic structures of order $`g`$ the edge $`(i,j)`$ belongs to, while $`s_{i,j}^{(g)}`$ is the number of possible cyclic structures of order $`g`$ that can be built given the degrees of the nodes.

We can now define, for every $`g`$, a detection algorithm that works exactly as the GN method with the difference that, at every step, the removed edges are those with the smallest value of $`\stackrel{~}{C}_{i,j}^{(g)}`$. By considering increasing values of $`g`$, one can smoothly interpolate between a local and a non local algorithm. Notice that the definition of $`\stackrel{~}{C}_{i,j}^{(g)}`$ guarantees that nodes with only one connection are not considered as isolated communities by the algorithm, since for their unique edge $`\stackrel{~}{C}_{i,j}^{(g)}`$ is infinite.

We have checked the accuracy of this new algorithm by comparing its performance with the GN method. Fig. 2 reports the results for the artificial test graph with four communities. It turns out that, with respect to the strong definition, the new algorithm is as accurate as the GN one, both in the case of cycles of order $`g=3`$ (triangles) and $`g=4`$ (squares).

On the other hand, for the weak definition the best accuracy is achieved with the new algorithm with $`g=4`$. Another test is performed by considering the examples of social networks already studied by GN. Fig. 3 shows the trees resulting from the application of the GN algorithm and of the $`g=4`$ algorithm to the network of college football teams. Again the results are very similar, indicating that the local algorithm captures well the presence of communities in that network.

Additional insight into the relationship between the GN algorithm and this new algorithm based on edge clustering is provided by Fig. 4, where the edge betweenness is plotted versus $`\stackrel{~}{C}_{i,j}^{(4)}`$ for each edge of the graph of scientific collaborations studied in Ref. Newman01bis . It is clear that an anticorrelation exists between the two quantities: edges with low values of $`\stackrel{~}{C}_{i,j}^{(4)}`$ tend to have high values of betweenness. The anticorrelation is not perfect: the edge with minimum $`\stackrel{~}{C}_{i,j}^{(4)}`$ is not the one with maximal betweenness. Therefore we expect the two algorithms to yield similar community structures though not perfectly coinciding.

In this framework it is important to recall the parallel drawn by Ronald S. Burt burt  between betweenness centrality and the so-called redundancy. The definition of redundancy is very close to that of node clustering and, much in the spirit of our work, Burt first pointed out that nodes that belong to few loops are central in the betweenness sense.

Let us now turn our attention to the question of the computational efficiency of our local algorithm. One can roughly estimate the scaling of the computational time as follows. When an edge is removed one has to check whether the whole system has been separated in disconnected components and update the value of $`\stackrel{~}{C}_{i,j}^{(g)}`$ in a small neighborhood of the removed edge. The first operation requires a time of the order of $`M`$, the total number of edges present in the network, while the time required by the second operation does not scale with $`M`$. Since this operation has to be repeated for all edges we can estimate the scaling of the total time as $`aM+bM^2`$. We thus expect computational time to be linearly dependent on $`M`$ for small systems and to cross over to an $`M^2`$ regime for large sizes.

We have measured the velocity of the new algorithm by computing the time needed to generate the whole tree for a random graph of increasing size $`N`$ and fixed average degree ($`MN`$). Results, reported in Fig. 5, confirm that the algorithm based on the computation of the edge clustering coefficient is much faster than the one based on edge betweenness, both for $`g=3`$ and $`g=4`$. The crossover between the initial linear dependence on $`N`$ to the $`N^2`$ growth for the local algorithm is evident for $`g=3`$.

## V Community structure in a network of scientific collaborations

In this section we consider an application of the fast algorithm discussed in the previous section to a network of scientific collaborations. In particular we have considered the network of scientists who signed at least one paper submitted to the E-print Archive relative to Condensed Matter in the period 1995-1999 (http://xxx.lanl.gov/archive/cond-mat). Data have been kindly provided by Mark Newman.

The network includes $`15616`$ nodes (scientists) but it has a giant connected component which includes only $`N=12722`$ scientists. We have focused our attention on this giant component and we have applied to it our algorithm for $`g=3`$ and $`g=4`$. The time needed for the generation of the dendrogram for $`g=3`$ is about 3 minutes on a desktop computer with 800 MHz CPU. The algorithm detects at the same time the communities satisfying the weak and the strong criteria. At the end of the procedure one obtains, for every value of $`g`$, a list of communities identified in a weak sense and another list of communities identified in a strong sense. By definition the second list is a subset of the first one. Figure 6 reports the size distributions (for $`g=3`$ and $`g=4`$) of all the communities identified in a weak sense. These distributions feature a power-law behaviour $`P(S)S^\tau `$ with $`\tau 2`$, an indication of the self-similar community structure of this network guido\_paolo . It is worth mentioning that the exponent we observe coincides with the one obtained, using the GN algorithm, in the framework of a recently proposed model of social network formation boguna .

For what concerns a more detailed analysis of the communities found, validation of the results is far from trivial, since there is no quantitative criterion to assess their accuracy. One may directly inspect the dendrogram in order to answer questions like: are the communities representative of real collaborations between the corresponding scientists? Do they identify specific research areas? Would a generic scientist agree about his or her belonging to a given community? Obviously all these questions cannot be answered in a definitive and quantitative way. We have partially followed this path and we have checked several subsets of the network at different levels in the hierarchy. To the best of our knowledge the results seem to us reasonable. Of course this does not represent a proof of the efficiency of the algorithm. We refer the reader to the detailed results of our analysis that we make available as additional supporting information 333The list of all the communities found with our algorithm for the scientific collaboration network is available on request (castella@pil.phys.uniroma1.it).

## VI Conclusions

The detection of the community structure in large complex networks is a promising field of research with many open challenges. The concept of community is qualitatively intuitive. However, to analyze a network it is necessary to specify quantitatively and unambiguously what a community is. Once a definition is given it is in principle possible to determine all subgraphs of a given network that fulfill the definition. However this task is in practice computationally out of reach even for small systems. Therefore the search for the community structure has generally a more limited goal: selecting, among all possible communities, a subset of them organized hierarchically, a dendrogram. Divisive and agglomerative algorithms carry out this task. A comparison of the performances of such algorithms is non trivial. In some simple cases, as the artificial graph with four subsets considered above, it is possible to assess quantitatively the validity of the results. In other cases, like the network of scientific collaborations, no quantitative measure exists to decide, given a precise definition of community, how good a dendrogram is. Typically, one may check whether the results appear sensible. However this is far from objective, being mediated by the observer’s own perception and by his/her intuitive concept of community.

In this work we have proposed two improvements in the construction of the dendrogram. First, we have devised a way to implement, in a generic divisive algorithm, a quantitative definition of community. In this way algorithms become fully self-contained, i.e. they do not need non-topological input to generate the dendrogram. Then we have introduced a new divisive algorithm, which is based on local quantities, and therefore extremely fast. Both these improvements have been tested successfully in controlled cases. The analysis of the large network of scientific collaborations gives results that appear reasonable. However, it is clear that, as discussed above, this statement is subjective and cannot be made at present more precise. Definitely, a quantitative measure for the evaluation of dendrograms would be a major step forward in this field.

At this point a remark is in order. So far we have only discussed examples of the so-called social networks. It has been shown newman-soc-nonsoc  that social networks substantially differ from other types of networks, namely technological or biological networks. Among other differences, they exhibit a positive correlation between the degree of adjacent vertices (assortativity) while most non-social networks are disassortative. While these results are consistent with our and other’s findings about community structures in social networks, they put into question the very existence of a community structure in non-social networks and the possibility of detecting it with the existing algorithms. From the perspective of our local algorithm, which relies on the existence of closed loops, disassortative networks could be in principle problematic, due to the small number of short cycles. However, interesting insight comes from the study of the loops of arbitrary order cal-pas-ves . In particular for four different types of networks (two of them social and assortative and two non-social and disassortative) measured values for the so-called average grid coefficient (the extension of the concept of clustering coefficient to cycles of order four) are two to four order of magnitude larger than the corresponding coefficients of a random graph with the same average degree and size $`N`$. This argues in favor of the presence of some sort of hierarchical structure and well-defined communities also in disassortative networks. It also hints that our algorithm could be fruitfully applied also to non-social (disassortative) networks, although future work will be needed in this direction.

We believe that the new elements presented in this paper can be of great help in the analysis of networks. On the one hand, the implementation of a quantitative definition of community makes algorithms self-contained and allows the analysis of the community structure based only on the network topology. On the other hand, the introduction of a new class of local and fast algorithms could open the way for applications to large-scale systems.

Acknowledgments: The authors are very grateful to Mark Newman for providing data on the networks of college football teams and of scientific collaborations. They also wish to thank Alain Barrat for useful suggestions and discussions.

Figure legends

* A simple network (left) and the corresponding dendrogram (right).
* Test of the efficiency of the different algorithms in the analysis of the artificial graph with four communities. The construction of the graph is described in text. Here $`N=128`$ and $`p_{in}`$ is changed with $`p_{out}`$ in order to keep the average degree equal to $`16`$. (Left) Strong definition: fraction of successes for the different algorithms compared with the analytical probability that four communities are actually defined. (Right) Weak definition: in addition to the same quantities plotted in the left graph, here we report, for every algorithm, the fraction $`f`$ of nodes not correctly classified.
* Plot of the dendrograms for the network of college football teams, obtained using the Girvan-Newman algorithm (left) and our new algorithm with $`g=4`$ (right). Different symbols denote teams belonging to different conferences. In both cases the observed communities perfectly correspond to the conferences, with the exception of the 6 members of the Independent conference, which are misclassified.
* Edge betweenness vs. the modified edge clustering coefficient $`\stackrel{~}{C}_{i,j}^{(4)}`$, for the network of scientific collaborations considered in section V. Each dot represents an edge in the network. See section V for details.
* Plot of the average time (in seconds) needed to analyze a random graph of $`N`$ nodes and fixed average degree $`k=5`$. The time refers to the construction of the full tree down to single nodes (Fig. 1). No criterion to validate the communities was imposed. The runs are performed on a desktop computer with a 800 MHz CPU.
* Normalized size distribution of all the communities of scientists identified in a weak sense by the algorithm described in section IV for $`g=3`$ (circles) and $`g=4`$ (squares). In both cases the behaviour is well reproduced by a power law with exponent $`2`$.


text:
# High Energy Dynamics of the Single Impurity Anderson Model

## Abstract

The quantitative control of the dynamic correlations of single impurity Anderson models is essential in several very active fields. We analyze the one-particle Green function with a constant energy resolution by dynamic density-matrix renormalization. In contrast to other approaches, sharp dominant resonances at high energies are found. Their origin and importance are discussed.

Single impurity models are at the very basis of the description of strong correlation phenomena. Landmarks are the Kondo problem kondo64  and the single impurity Anderson model (SIAM) ander61 , for a review see Ref. hewso93, .

The interest in the quantitative analysis of SIAMs has been intensified considerably by the advent of a systematic mapping of strongly correlated lattice models to effective SIAMs supplemented by a self-consistency condition. This is the key point of dynamic mean-field theory (DMFT) jarre92 ; georg92a  which is based on an appropriate scaling of the non-local parts of the Hamiltonian metzn89a ; mulle89a , for reviews see Ref. prusc95, ; georg96, . In recent years, the DMFT is applied very successfully in combination with ab-initio density functional calculations anisi97 ; licht98 . In this way, the unbiased knowledge about the bands could be enhanced by the inclusion of interaction effects between the excited quasi-particles. It turned out that the combination of density functional results and DMFT makes the quantitative understanding of spectroscopic data possible held01b .

So far, the methods applied to the SIAM were designed to capture the low-energy physics, in particular the fixed points and the thermodynamics krish80a ; schlo82 . The numerical renormalization group (NRG) was later extended to calculate also dynamic, i.e., spectral information. It provides reliable data on the scale of the Kondo temperatures $`T_\mathrm{K}`$ , see Ref. hewso93, ; bulla00a,  and references therein. On larger scales, the energy resolution is less well-controlled.

But in various applications the behavior at higher energies is important to achieve quantitative accuracy. For instance, the self-consistency cycle of the DMFT mixes modes at all energies. Hence, excellent quantitative control over the dynamics at high energies is indispensable, even if finally only the behavior at low energies matters.

Another application is the optical control of isolated $`S=1/2`$ impurities or quantum dots coupled to narrow bands. If the impurities differ so that the energy between the singly occupied ground state and the excited double occupancy differs, they can be switched selectively from the ground state to the double occupancy (and back) by shining light at the resonant frequency onto the sample. The life-time of the double occupancy, i.e., the inverse line width of the resonance, determines how well the resonance condition has to be met, how selective the switching can be, and how stable the excited state is.

In view of the above, we perform a numerical investigation which aims to describe both the low-energy dynamics and the high-energy dynamics quantitatively. To this end, we use an energy resolution which is constant for all energies. Features at low energies are not as delicately resolved as by NRG, but in return features at high energies are much better under control. We apply the dynamic density-matrix renormalization (D-DMRG) hallb95 ; kuhne99a ; hovel00  to compute the one-particle propagator. The DMRG is a real-space approach white92 ; white93 ; pesch98  which works best for open boundary conditions so that it is particularly well-suited to treat impurity problems.

The model investigated at zero temperature is the symmetric Anderson model

$`={\displaystyle \underset{\sigma }{}}ϵ_dn_{d,\sigma }+Un_{d,}n_{d,}+V{\displaystyle \underset{\sigma }{}}(d_\sigma ^{}c_{0,\sigma }^{}+\text{h.c.})`$
$`+{\displaystyle \underset{n,\sigma }{}}\gamma _{n+1}\left(c_{n,\sigma }^{}c_{n+1,\sigma }^{}+\text{h.c.}\right)+{\displaystyle \underset{n,\sigma }{}}ϵ_nc_{n,\sigma }^{}c_{n,\sigma }^{}`$ (1)
with arbitrary density of states (DOS) $`\rho _0(\omega )`$ of the $`U=0`$ one-particle Green function $`G_0(\omega )`$ of the $`d`$-electron. The parametrization in (1) is chosen such that the coefficients $`(ϵ_n,\gamma _n)`$ are the continued fraction coefficients of the hybridization function

$$\mathrm{\Gamma }(\omega )=\frac{V^2}{\omega ϵ_0{\displaystyle \frac{\gamma _1^2}{\omega ϵ_1{\displaystyle \frac{\gamma _2^2}{\omega \mathrm{}}}}}}.$$
(2)
This model has particle-hole symmetry iff $`ϵ_d=U/2`$ and $`ϵ_n=0`$ for all $`N0`$. The representation as continued fraction viswa94  (see Fig. 1a) is optimum for the DMRG which is designed for chains. We look at a generic situation with finite band width $`W=2D`$. For simplicity we choose a $`\mathrm{\Gamma }(\omega )`$ with semi-elliptic DOS, i.e., $`\gamma _n=D/2`$. For $`V=D/2`$, the free DOS $`\rho _0(\omega )=2\sqrt{D^2\omega ^2}/(\pi D^2)`$ is also semi-elliptic.

The problem illustrated in Fig. 1a is mapped by two standard Jordan-Wigner transformations to two $`XY`$ spin $`1/2`$ chains, the $`S`$-chain and the $`T`$-chain. The $`S`$-chain results from the $``$ fermions, the $`T`$-chain from the $``$ fermions. They are coupled at site zero where the density-density coupling is mapped to the product of $`z`$-components. The resulting chain is depicted for the symmetric SIAM in Fig. 1b. The couplings are given by $`t_0=t_0^{}=V`$ and $`t_i=t_i^{}=\gamma _i`$ for $`i1`$. The mapping from fermions to spins avoids the fermionic Fock space which would imply numerically difficult long-range effects. The mapping makes the Hilbert space the direct product of the local Hilbert spaces at each site.

The DMRG can easily determine the ground state $`|0`$ and its energy $`E_0`$ for a finite chain. So the chain in Fig. 1b is truncated such that there are $`L`$ spins in the upper and in the lower part of the chain corresponding originally to a truncated bath of $`L1`$ fermions plus the impurity. The dynamic quantity we are interested in is the retarded Green function at zero temperature

$$G^>(\omega +i\eta )=0|𝐒_0^{}(\omega +i\eta (E_0))^1𝐒_0^+|0$$
(3)
where the superscript <sup>\></sup> implies that (3) represents only the part of the usual Green function at non-negative frequencies. In the symmetric case, the complete function is recovered by $`G(\omega +i\eta )=G^>(\omega +i\eta )G^>(\omega i\eta )`$. In the asymmetric case, $`G^<(\omega +i\eta )=0|𝐒_0^+(\omega +i\eta (E_0))^1𝐒_0^{}|0`$ must be determined separately whereby $`G(\omega +i\eta )=G^>(\omega +i\eta )G^<(\omega i\eta )`$ is obtained. We stress that $`G(\omega +i\eta )`$ is the fermionic propagator even though it is computed in terms of spins after the Jordan-Wigner mapping.

The key idea of the dynamic DMRG is to include the real and the imaginary part of the correction vector $`|c`$ in the target states of a standard DMRG algorithm kuhne99a ; hovel00 . The natural choice is $`|c=(\omega +i\eta (E_0))^1𝐒_0^+|0`$. The computation of $`|c`$ is numerically the most demanding step due to the inversion of an almost singular non-hermitean matrix. We prefer to stabilize this inversion by optimized algorithms freun93a  instead of using the variational approach proposed by Jeckelmann jecke02  which requires a minimization in a high-dimensional Hilbert space.

The numerical calculations cannot be performed for $`\eta =0`$. Even small values of $`\eta `$ are very time consuming. So we compute first $`G(\omega +i\eta )`$ at finite $`\eta `$. The spectral density $`\rho ^{(\eta )}(\omega )=\frac{1}{\pi }\mathrm{}G(\omega +i\eta )`$ can be seen as the actual spectral density $`\rho (\omega )`$ convoluted by the Lorentzian $`\rho _\text{L}(\omega )=(\eta /\pi )/(\omega ^2+\eta ^2)`$ of width $`\eta `$. Hence it is possible to retrieve $`\rho (\omega )`$ by deconvolution. A standard technique for deconvolution is Fourier transformation, realized best by fast Fourier transforms, division by $`\mathrm{exp}(\eta \tau )`$ plus low-pass filtering, and the inverse transform. A flexible alternative with similar properties is the explicit matrix inversion of the convolution procedure gebha03 .

In Fig. 2a, generic broadened spectral densities are plotted as they are computed by D-DMRG. Obviously, the value $`\rho ^{(\eta )}(0)`$ is not independent of $`U`$. Increasing the chain length $`L`$ does not lead to any significant change in the data (not shown). Fig. 2b displays the deconvoluted data. The deconvolution works very well except for some slight overshooting in regions where the spectral density changes rather abruptly. In particular, the value $`\rho (0)`$ is pinned to $`D/(2\pi V^2)`$ independent of $`U`$ as required by Friedel’s sum rule and the density of states rule lutti60b ; lutti61 ; ander91 ; hewso93 . We take this fact as convincing evidence for the reliability of the numerical algorithm.

The central peak at $`\omega =0`$ is the Abrikosov-Suhl resonance (ASR). For larger $`U`$ (smaller $`V`$) its width decreases rapidly so that the ASR is very difficult to resolve nishi03  unless more elaborate deconvolution schemes are devised raas04a . So a quantitative analysis of the ASR is postponed to future work.

For comparison, the thin dashed lines in Fig. 2b depict standard NRG data bulla98 ; bulla00a . For small frequencies NRG is well-controlled. Indeed, for $`|\omega |D/3`$, NRG and D-DMRG data agree excellently lending further support to the D-DMRG approach. Outside the band, the NRG spectra appear to be too wide due to the chosen constant broadening on a logarithmic mesh. This broadening does not account for the absence of states outside the bare band. The NRG does not possess intrinsic information about the peak widths. The position of the high energy peak in the raw NRG data, however, coincides with the D-DMRG result.

An increase in $`U`$ leads to the formation of Hubbard satellites below and above the free band (Fig. 2). They are situated at energies $`\omega _{\mathrm{up}/\mathrm{low}}=\pm (U/2+\delta _{\mathrm{shift}}),\delta _{\mathrm{shift}}>0`$ and become more pronounced on increasing $`U`$ in two ways. They capture more weight and they become sharper. For $`T_\text{K}0`$ the weight reaches 1/2, see Ref. ander91, . The sharpening has not been discussed quantitatively before although the extended non-crossing approximation prusc89  provides sharp satellites if they lie outside the bare bands, see e.g. Fig. 1 in Ref. ander91, . Recently, indications have occured gebha03  that other standard algorithms overestimate the width of the Hubbard satellites. The exaggerated width of the NRG data at high energies results from the Gaussian broadening of the order of the energy range sakai89 .

To investigate the line shapes of the satellites we plot them for various values of $`U`$ and $`V`$ in Fig. 3. The ASR at $`\omega =0`$ is not displayed since it is too much smeared out at $`\eta =0.1D`$ for larger values of the interaction. The shifts $`\delta _{\mathrm{shift}}`$ increase on increasing $`V`$; they decrease on growing interaction $`U`$. The widths behave qualitatively similar. A complete deconvolution suffers unfortunately from severe overshooting due to the sharpness of the resonance. To make the analysis nonetheless quantitative we fit the broadened data by Lorentzians plus an offset (Fig. 3). These fits work very well for large values of $`U`$ and not too large values of $`V`$.

To deduce the true width of the Hubbard satellite we assume that it is well described by a Lorentzian. The width $`\eta _{\mathrm{eff}}`$ of the convolution of two Lorentzians of widths $`\eta _1`$ and $`\eta _2`$ is $`\eta _{\mathrm{eff}}=\eta _1+\eta _2`$. From the effective widths $`\eta _{\mathrm{eff}}`$ we deduce the true half-widths at half-maximum (HWHM) of the Hubbard satellite by subtracting the artificial broadening $`\eta `$, i.e., HWHM $`=\eta _{\mathrm{eff}}\eta `$. In Fig. 4, the widths are depicted as function of $`U`$ and of $`V`$. The results show that the HWHM are proportional to $`V^4/U^2`$. The deviations for smaller widths must be attributed to the numerical constraints, e.g., finite $`\eta `$ and finite chain length $`L`$. The deviations for larger widths, mainly for larger values of $`V`$ and smaller values of $`U`$ result from the vicinity of the bare bands.

Fig. 5 displays the analogous analysis for the shifts $`\delta _{\mathrm{shift}}`$ of the satellite positions. Again, strong evidence for power law behavior is found, namely $`\delta _{\mathrm{shift}}V^2/U`$.

How can the above findings be understood? Let us start by the positions. The energy levels of isolated impurities, i.e., $`V=0`$ are at $`\pm U/2`$ schri66 . Switching on $`V`$ mixes the impurity levels with the bath states which lie in the interval $`[D,D]`$. If $`U`$ is large compared to $`D`$ second order perturbation in $`V`$ implies that the impurity levels are repelled from the bath states. The shift should be of the order of $`J=8V^2/U`$, see Eq. (11) in Ref. schri66, , which agrees nicely with the power laws in Fig. 5.

The widths of the satellites have been discussed quantitatively when they lie within the bare band logan98 . If the satellites lie outside, perturbation theory in $`U`$, to second order or the random phase approximation, implies that a finite width is to be expected at least for $`U<6D`$. But the reasoning in powers of $`U`$ is not helpful for $`U>6D`$ and it does not explain the power laws found. So we return to perturbing in powers of $`V`$. The impurity levels mix with particle-hole excitations in the bands, see Eq. (10a) in Ref. schri66, . In the symmetric case the doubly occupied electron and hole state are degenerate so that mixing with particle-particle (or hole-hole) states matters also, see Eq. (12) in Ref. schri66, . The mixing is of order $`J=8V^2/U`$. So Fermi’s golden rule implies a life time of $`J^2N_0`$ where $`N_0`$ measures the density of states with which the impurity level mixes; $`N_0`$ is of the order of $`D^1`$. Indeed, HWHM $`J^2`$ explains conclusively the data of Fig. 4.

So far, the width of the Hubbard satellites for $`U>2D`$ was extracted under the assumption that the satellites are Lorentzians. Further investigations of the line shape are urgently called for. Numerically, improvements of the resolution are necessary to determine the line shape of the satellites explicitly. Analytically, the quantitative argument for the widths must be supplemented by an explicit calculation of the line shapes for $`V^2/U0`$.

In summary, we have investigated the dynamic propagator of the SIAM by D-DMRG. This powerful large-scale algorithm provides information with a constant energy resolution. Up to moderate interactions $`U2D`$, deconvolution yields the explicit spectral densities. For larger interactions, the width of sharp resonances can be extracted by fitting Lorentzians. In particular, we analyzed the positions and widths of the Hubbard satellites. The shifts are of order $`V^2/U`$ due to level repulsion; the line widths are of order $`V^4/U^2`$.

Especially the sharpness of Hubbard peaks is missed by other zero temperature algorithms for the SIAM. Hence the D-DMRG is a very valuable complementary tool. Position and width of the Hubbard satellites are important for several widely used applications, e.g., in the self-consistency cycle of the DMFT.

Acknowledgements We thank R. Bulla, L. Craco, M. Grüninger, T. Hövelborn, D. Logan, H. Monien, and E. Müller-Hartmann for helpful discussions and the DFG for financial support (Uh 90/3-1 and SFB 608).


text:
# 1 Introduction

## 1 Introduction

The production of heavy hadrons ($`H`$) in $`e^+e^{}`$annihilation provides a laboratory for the study of heavy-quark ($`Q`$) jet fragmentation. This is commonly characterized in terms of the observable $`x_H`$ $``$ $`2E_H/\sqrt{s}`$, where $`E_H`$ is the energy of a $`B`$ or $`D`$ hadron containing a $`b`$ or $`c`$ quark, respectively, and $`\sqrt{s}`$ is the c.m. energy. The distribution of $`x_H`$, $`D(x_H)`$, is conventionally referred to as the heavy-quark ‘fragmentation function’***Unless stated otherwise, in these studies we do not distinguish between hadrons and antihadrons..

In recent publications we presented  the results of a new method for reconstructing $`B`$-hadron decays, and the $`B`$ energy, inclusively, using only charged tracks, in the SLD experiment at SLAC. We used the upgraded charge-coupled device (CCD) vertex detector, installed in 1996, to reconstruct $`B`$-decay vertices in $`Z^0`$ decays with high efficiency and purity. Combined with the micron-sized SLC interaction point (IP), our precise vertexing allowed us to reconstruct the total transverse momentum of the tracks from $`B`$-decays, and therefore the transverse momentum and mass associated with the neutral particles in the $`B`$-decays. This allowed us to reconstruct accurately the energy of $`B`$ hadrons. These studies yielded the most precise measurement of the $`b`$-quark fragmentation function, and allowed us to test models of heavy-quark fragmentation. Of the 9 models tested, only 4 were consistent with our precision data at better than the 1% level based on a $`\chi ^2`$ probability. This allowed us to reduce the model-dependent systematic uncertainty on the $`b`$-quark fragmentation function.

We have extended these studies and applied similar ‘topological’ vertexing techniques to tag events in which we reconstructed the energies of both leading $`B`$ hadrons produced via $`e^+e^{}`$$``$$`b\overline{b}`$$``$$`B\overline{B}+X`$. We measured the moments of the single-inclusive $`B`$-hadron scaled-energy distribution $`dN/dx_B`$:

$$D_ix_B^{i1}\frac{1}{N_s}\frac{dN_s}{dx_B}𝑑x_B$$
(1)
as well as the moments of the double-inclusive scaled-energy distribution:

$$D_{ij}(\varphi )x_{B1}^{i1}x_{B2}^{j1}\frac{1}{N_d}\frac{d^3N_d}{dx_{B1}dx_{B2}d\mathrm{cos}\varphi }𝑑x_{B1}𝑑x_{B2},$$
(2)
where $`x_{B1}`$ and $`x_{B2}`$ are the scaled energies of the two $`B`$ hadrons and the label is arbitrary, $`\varphi `$ is the angle between their flight directions, and $`i`$ and $`j`$ are integers $``$1. We formed the normalised moments:

$$G_{ij}(\varphi )D_{ij}(\varphi )/(D_iD_j)$$
(3)
Following the method proposed in  we used these quantities to test the ansatz of factorisation as applied to perturbative Quantum Chromodynamics (pQCD) calculations of $`e^+e^{}`$$``$$`b\overline{b}`$events.

Due to ‘soft’, or non-perturbative, effects, $`D_i`$ and $`D_{ij}`$ cannot be predicted absolutely in pQCD. Rather, the respective pQCD calculation must be folded with models of the non-perturbative (‘hadronisation’) process in order to derive predictions that can be compared with experimental data. However, provided that the ansatz of factorisation holds , namely that calculation of the perturbative and non-perturbative phases can be separated by (an arbitrary) factorisation scale $`\mu _F`$, the dependence on $`\mu _F`$ cancels in Eq. 3 and hence $`G_{ij}`$ can be calculated absolutely in pQCD, up to possible ‘higher twist’ effects of order $`1/\sqrt{s}`$, with no dependence on hadronisation models. Comparison of the measured $`G_{ij}`$ with pQCD predictions hence allows both a test of this ansatz and of the perturbative calculations. The $`G_{ij}`$ can be derived at leading order (LO) in pQCD using the numerical results in . In addition, next-to-leading order (NLO) predictions for $`G_{ij}`$ have been calculated recently .

In Section 2 we describe the detector and the selection of $`e^+e^{}`$$``$hadrons events used in this analysis. We present in Section 3 the first measurement of the double-inclusive $`B`$-hadron energy distribution and, in Section 4, of the normalised moments Eq. 3. In Section 5 we describe the estimation of the errors on our measurements. Finally, in Section 6, we test the ansatz of factorisation and compare our data with the pQCD predictions.

## 2 Apparatus and Hadronic Event Selection

This analysis is based on roughly 400,000 hadronic $`Z^0`$ events produced in e+e- annihilations at a mean center-of-mass energy of $`\sqrt{s}=91.28`$ GeV at the SLAC Linear Collider (SLC), and recorded in the SLC Large Detector (SLD) between 1996 and 1998. A general description of the SLD can be found elsewhere . The trigger and initial selection criteria for hadronic $`Z^0`$ decays are described in Ref. . This analysis used charged tracks measured in the Central Drift Chamber (CDC)  and in the upgraded Vertex Detector (VXD3) . Momentum measurement was enabled by a uniform axial magnetic field of 0.6T. The CDC and VXD3 gave a momentum resolution of $`\sigma _p_{}/p_{}`$ = $`0.010.0026p_{}`$, where $`p_{}`$ is the track momentum transverse to the beam axis in GeV/$`c`$. In the plane normal to the beamline the centroid of the micron-sized SLC IP was reconstructed from tracks in sets of approximately thirty sequential hadronic $`Z^0`$ decays to a precision of $`\sigma _{IP}^{r\varphi }4`$ $`\mu `$m. The IP position along the beam axis was determined event by event, using charged tracks, with a resolution of $`\sigma _{IP}^z`$ $``$ 20 $`\mu `$m. Including the uncertainty on the IP position, the resolution on the charged-track impact parameter ($`d`$) projected in the plane perpendicular to the beamline was $`\sigma _d^{r\varphi }`$ = 8$``$33/$`(p\mathrm{sin}^{3/2}\theta )`$ $`\mu `$m, and the resolution in the plane containing the beam axis was $`\sigma _d^z`$ = 10$``$33/$`(p\mathrm{sin}^{3/2}\theta )`$ $`\mu `$m, where $`\theta `$ is the track polar angle with respect to the beamline.

A set of cuts was applied to the data to select well-measured tracks and events well contained within the detector acceptance. Charged tracks were required to have a distance of closest approach transverse to the beam axis within 5 cm, and within 10 cm along the axis from the measured IP, as well as $`|\mathrm{cos}\theta |<0.87`$, and $`p_{}>0.15`$ GeV/c. Events were required to have a minimum of five such tracks, a thrust axis  polar angle w.r.t. the beamline, $`\theta _T`$, within $`|\mathrm{cos}\theta _T|<0.80`$, and a charged visible energy $`E_{vis}`$ of at least 20 GeV, which was calculated from the selected tracks assigned the charged pion mass. The efficiency for selecting a well-contained $`Z^0q\overline{q}(g)`$ event was estimated to be above 97% independent of quark flavor. The selected sample comprised 313,447 events, with an estimated $`0.10\pm 0.05\%`$ background contribution dominated by $`Z^0\tau ^+\tau ^{}`$ events.

For the purpose of estimating the efficiency and purity of the selection procedures we made use of a detailed Monte Carlo (MC) simulation of the detector. The JETSET 7.4  event generator was used, with parameter values tuned to hadronic $`e^+e^{}`$annihilation data , combined with a simulation of $`B`$-hadron decays tuned  to $`\mathrm{{\rm Y}}(4S)`$ data and a simulation of the SLD based on GEANT 3.21 . Inclusive distributions of single-particle and event-topology observables in hadronic events were found to be well described by the simulation . Uncertainties in the simulation were taken into account in the systematic errors (Section 5).

## 3 $`B`$-Hadron Selection and Energy Measurement

The event sample for this analysis was selected using a ‘topological’ vertexing technique based on the detection and measurement of charged tracks, which is described in detail in Refs. . We considered events in which we found decay vertices corresponding to both the leading $`B`$ and $`\overline{B}`$ hadrons. First, the Durham algorithm  was applied to the selected hadronic events, with a $`y_c`$ parameter value of 0.015, in order to define a jet structure in each event. We found that this algorithm and $`y_c`$ value minimized the number of $`B`$ (and $`D`$) decay tracks assigned to the wrong jet. This is an important feature for our analysis because we used vertex-related variables derived only from charged tracks. Events containing 2, 3, or 4 jets were retained for further analysis.

In each selected event, the vertexing algorithm was applied to the set of tracks in each jet. Vertices consistent with photon conversions or $`K^0`$ or $`\mathrm{\Lambda }^0`$ decays were discarded. Events were retained in which a vertex was found in exactly two jets. 35,137 events were selected, of which 89.4% were estimated to be of $`b\overline{b}`$origin. The efficiency for selecting true $`b\overline{b}`$events was estimated to 36.3%.

The large masses of the $`B`$ hadrons relative to light-flavor hadrons make it possible to distinguish $`B`$-hadron decay vertices from those vertices found in events of lighter flavors using the vertex invariant mass, $`M`$. However, due to those particles missed from the vertex, which are mainly neutrals, $`M`$ cannot be fully determined. $`M`$ can be written

$$M=\sqrt{M_{ch}^2+P_t^2+P_{chl}^2}+\sqrt{M_0^2+P_t^2+P_{0l}^2}$$
(4)
where $`M_{ch}`$ and $`M_0`$ are the total invariant masses of the set of vertex-associated tracks and the set of missing particles, respectively. $`P_t`$ is the total charged-track momentum transverse to the $`B`$ flight direction, which, by momentum conservation, is identical to the transverse momentum of the set of missing particles. $`P_{chl}`$ and $`P_{0l}`$ are the respective momenta along the $`B`$ flight direction, which we take to be the vector joining the IP to the vertex.

The lower bound for the mass of the decaying hadron, the ‘$`P_t`$-corrected vertex mass’ ,

$$M_{Pt}=\sqrt{M_{ch}^2+P_t^2}+|P_t|$$
(5)
was used as the variable for selecting $`B`$ hadrons. Figure 1 shows the distribution of $`M_{Pt}`$ for vertices in the selected event sample, and the corresponding simulated distribution. Events were selected that contained at least one vertex with $`M_{Pt}`$ $`>`$ 2.0 GeV/$`c^2`$ and $`M_{Pt}2\times M_{ch}`$. The latter cut was found to reduce the contamination from fake vertices in light-quark events.

In order to improve the $`b\overline{b}`$purity of the sample, events were selected in which both vertices had a flight length, $`d_{vtx}`$, such that $`0.1<d_{vtx}<2.3`$ cm; in which at least one vertex contained two ‘significant’ tracks, i.e. tracks with a normalised impact-parameter significance, $`d/\sigma _d`$, of at least 2 units; and in which the angle between the vertex flight vectors, $`\varphi `$, satisfied cos$`\varphi `$ $`<`$ 0.99. The last cut was effective at removing events in which either a gluon had split into heavy quarks, or the jet-finder had artificially split a single heavy-quark jet into two jets.

The energy of each $`B`$ hadron, $`E_B`$, can be expressed as the sum of the reconstructed-vertex energy, $`E_{ch}`$, and the energy of those particles not associated with the vertex, $`E_0`$. We can write

$$E_0^2=M_0^2+P_t^2+P_{0l}^2$$
(6)
The two unknowns, $`M_0`$ and $`P_{0l}`$, must be found in order to obtain $`E_0`$. One kinematic constraint can be obtained by imposing the $`B`$-hadron mass on the vertex, $`M_B^2=E_B^2P_B^2`$, where $`P_B=P_{chl}+P_{0l}`$ is the total momentum of the $`B`$ hadron. From Equation (4) we derive the following inequality,

$$\sqrt{M_{ch}^2+P_t^2}+\sqrt{M_0^2+P_t^2}M_B,$$
(7)
where equality holds in the limit where both $`P_{0l}`$ and $`P_{chl}`$ vanish in the $`B`$ hadron rest frame. Equation (7) effectively sets an upper bound on $`M_0`$, and a lower bound is given by zero:

$$0M_0^2M_{0max}^2,$$
(8)
where

$$M_{0max}^2=M_B^22M_B\sqrt{M_{ch}^2+P_t^2}+M_{ch}^2.$$
(9)
Because $`M_0`$ peaks near $`M_{0max}`$ we set $`M_0^2`$ = $`M_{0max}^2`$ if $`M_{0max}^2`$ $``$0, and $`M_0^2`$ = 0 if $`M_{0max}^2`$ $`<`$0. We calculated $`P_{0l}`$:

$$P_{0l}=\frac{M_B^2(M_{ch}^2+P_t^2)(M_0^2+P_t^2)}{2(M_{ch}^2+P_t^2)}P_{chl},$$
(10)
and hence $`E_0`$ (Equation (6)). We then reconstructed the $`B`$ hadron energy, $`E_B^{rec}=E_0+E_{ch}`$. Events were retained in which both reconstructed $`B`$ energies satisfied $`E_B^{rec}<60`$ GeV. A final sample of 19,809 events was obtained with an estimated $`b\overline{b}`$selection efficiency of 21.7% and a background contribution of only 0.23%, which was almost entirely from $`Z^0`$ $``$$`c\overline{c}`$events.

The energy resolution of the final $`B`$ sample is shown in Fig. 2, where we plot the normalised residual on $`E_B`$: $`(E_B^{true}E_B^{rec})/E_B^{true}`$. The distribution was fitted to a double Gaussian function for which the mean positions, widths and normalisations were allowed to vary. 79.5% of the population lies in the ‘core’ Gaussian, of width 21.3%; the remaining population is characterised by a ‘tail’ Gaussian of width 31.3%.

The angular resolution was similarly investigated and is shown in Fig. 3, where we plot the residual on $`\varphi `$: $`\varphi ^{true}\varphi ^{rec}`$. The vast majority of angles are reconstructed to better than 5 mrad. The tail in the resolution function corresponds to $`B`$ decays with shorter decay lengths, where the vertex position error had a larger relative effect on the determination of $`\varphi `$.

The double-inclusive distribution of raw $`B`$-hadron energies is shown in Figure 4. Since we do not distinguish between $`B`$ and $`\overline{B}`$ hadrons, the reconstructed $`B`$-decay vertices were labelled arbitrarily ‘1’ and ‘2’ on an event-by-event basis.

We divided $`E_B^{rec}`$ by the beam energy, $`E_{beam}=\sqrt{s}/2`$, to obtain the reconstructed scaled $`B`$-hadron energy, $`x_B^{rec}=E_B^{rec}/E_{beam}`$.

## 4 Angle-Dependent $`B`$-$`\overline{B}`$ Energy Moments

In each event we quantified the correlations between the two $`B`$ hadrons in terms of the angle dependent scaled-energy moments proposed in . We formed the single-inclusive $`B`$-energy distribution and evaluated the moments (Eq. 1) from the raw measured distribution:

$$D_i^{rec}=\frac{1}{2N}\mathrm{\Sigma }_k(x_B^{rec})^{i1}$$
(11)
where the index $`k`$ ($`1k20`$) runs over the bins of $`x_B^{rec}`$ and $`N`$ is the number of events in the sample. Similarly, from the double-inclusive $`B`$-energy distribution (Fig. 4) we evaluated in each cos$`\varphi `$ bin the double moments (Eq. 2):

$$D_{ij}^{rec}(\varphi )=\frac{1}{N}\mathrm{\Sigma }_{N(\varphi )}(x_{B1}^{rec})^{i1}(x_{B2}^{rec})^{j1}$$
(12)
where the sum extends over the set of events in each cos$`\varphi `$ bin. The normalised moments $`G_{ij}`$ (Eq. 3) were evaluated:

$$G_{ij}^{rec}(\varphi )=\frac{D_{ij}^{rec}(\varphi )}{D_i^{rec}D_j^{rec}}$$
(13)
The first six moments, $`i`$ = 1,2,3 and $`j`$ = 1,$`\mathrm{}`$,$`i`$ are shown in Fig. 5. The bin centers were defined by taking the average value of cos$`\varphi `$ within each bin.

Also shown in Fig. 5 is a comparison with the simulated normalised moments; the simulation reproduces the data. Given that we showed previously  that the Peterson function implemented in our simulation does not provide a good description of the $`b`$-quark fragmentation function, this agreement may naively appear to be surprising. However, if, as proposed in , the non-perturbative contributions to the normalised quantity $`G_{ij}`$ cancel, the agreement should be excellent, as observed.

We used our simulated event sample to correct for the effects of the detector acceptance, the efficiency of the technique for reconstructing $`B`$-hadron decays, the energy resolution, and bin migrations caused by the finite angular resolution. We defined a binwise correction factor:

$$F_{ij}^{MC}(\varphi )\frac{G_{ij}^{gen}(\varphi )|_{MC}}{G_{ij}^{rec}(\varphi )|_{MC}}$$
(14)
where $`G_{ij}^{gen}(\varphi )|_{MC}`$ is the normalised moment calculated using generated true $`e^+e^{}`$$``$$`B\overline{B}`$ \+ $`X`$ events, and $`G_{ij}^{rec}(\varphi )|_{MC}`$ is the corresponding moment calculated after simulation of the detector and application of the same analysis as applied to the data. For cos$`\varphi `$ $``$ $`1`$ the value of $`F_{ij}`$ is close to unity. As cos$`\varphi `$ increases towards 1, $`F_{ij}`$ rises monotonically to approximately 8 ($`F_{11}`$) or 2 ($`F_{33}`$). The increase with cos$`\varphi `$ reflects our decreasing efficiency to select $`B\overline{B}`$ events as the angle between the two $`B`$-decay vertices becomes smaller and the $`B`$ energies decrease. For a given cos$`\varphi `$ bin, $`F_{ij}`$ decreases as $`i`$ and $`j`$ increase due to the effective weighting of this efficiency by $`x_B^{i1}x_B^{j1}`$.

We derived the true normalised moments:

$$G_{ij}(\varphi )=F_{ij}^{MC}(\varphi )G_{ij}^{rec}(\varphi )$$
(15)
which are shown in Fig 6. The uncertainties associated with this correction procedure are discussed in the next section.

## 5 Error Estimation

We used our simulation to estimate the statistical error on the moments. The simulated event sample was divided into subsamples, each comprising the same number of events as in the data sample. The entire analysis was performed on each of these subsamples. Within each cos$`\varphi `$ bin the r.m.s. deviation of the ensemble of results was evaluated and taken as the statistical error in that bin. For a given pair of $`i`$ and $`j`$ values the errors on the particular moment are not correlated between different bins of cos$`\varphi `$, but within a given cos$`\varphi `$ bin the errors on the different moments are correlated, since the same set of events was used.

We considered sources of systematic uncertainty that potentially affect our measurements. These may be divided into uncertainties in modeling the detector and uncertainties on experimental measurements serving as input parameters to the underlying physics modeling. For these studies our simulation was used.

Since our energy reconstruction technique is strongly dependent on charged-track properties, four sources of uncertainty were investigated: our simulated tracking efficiency, transverse momentum resolution, and the resolutions on the track impact parameter and polar angles. In each case the simulation was corrected so as to reproduce the data. The full analysis was repeated on the data, and half the difference between the results obtained using the corrected and uncorrected simulations was taken as a symmetric systematic error.

A large number of measured quantities relating to the production and decay of charm and bottom hadrons are used as input to our simulation. In $`b\overline{b}`$events we considered the uncertainties on: the branching fraction for $`Z^0`$ $``$$`b\overline{b}`$; the rates of production of $`B^\pm `$, $`B^0`$ and $`B_s^0`$ mesons, and $`B`$ baryons; the lifetimes of $`B`$ mesons and baryons; and the average $`B`$-hadron decay charged multiplicity. In $`c\overline{c}`$events we considered the uncertainties on: the branching fraction for $`Z^0`$ $``$$`c\overline{c}`$; the charmed hadron lifetimes, the charged multiplicity of charmed hadron decays, the production of $`K^0`$ from charmed hadron decays, and the fraction of charmed hadron decays containing no $`\pi ^0`$s. We also considered the rates of production of secondary $`b\overline{b}`$and $`c\overline{c}`$from gluon splitting. The uncertainty on the world-average value  of each quantity was used to rederive the corrected $`G_{ij}(\varphi )`$. In each bin of $`G_{ij}(\varphi )`$ the deviation between the rederived and standard values was taken as an estimate of the corresponding systematic error. Other relevant systematic effects such as variation of the event selection cuts and the assumed $`B`$-hadron mass were found to be very small.

The error associated with the choice of $`b`$-quark fragmentation function used in the simulation was estimated by rederiving the corrected data using in turn the UCLA, Kartvelishvili and Bowler fragmentation functions . In each cos$`\varphi `$ bin the r.m.s. deviation of the results w.r.t. the standard value, using the Peterson function, was taken as an estimate of the systematic error. This error was typically much smaller than that arising from the other error sources.

For each systematic error investigated, the shape of the angular dependence of the moments was not significantly changed, only the overall normalization was affected. For each moment, in each cos$`\varphi `$ bin the dominant errors were typically those related to the uncertainties on the charged-track properties. The errors due to charm and bottom physics modeling were typically an order of magnitude smaller. In each cos$`\varphi `$ bin all sources of systematic uncertainty were added in quadrature to obtain the total systematic error. The fully-corrected $`G_{ij}`$, with statistical and systematic errors added in quadrature, are shown in Fig. 6.

## 6 Comparison with Perturbative QCD Calculations

The fully-corrected data were compared (Fig. 6) with a recent calculation  of the normalised moments complete at NLO in pQCD. Both the LO and full NLO calculations are shown in Fig. 6; the calculations assume an $`\alpha _s(M_Z^2)`$value of 0.120 and a pole $`b`$-quark mass value of 5.0 GeV/$`c^2`$. It can be seen that the difference between the two is relatively small, and that the LO calculation lies systematically slightly below the NLO calculation. The small size of the NLO relative to the LO contributions is an indication that the normalised moments are perturbatively robust observables. For each moment shown, the LO calculation undershoots the data. The NLO calculation reproduces the data across the full range of cos$`\varphi `$.

This comparison does not rely on any convolution of the pQCD calculations with models of the non-perturbative hadronisation process. Hence the excellent agreement between the pQCD calculations and the data verifies the ansatz of factorization between the perturbative and non-perturbative phases that forms the basis for the pQCD calculation of heavy-hadron properties.

## Acknowledgements

We thank the personnel of the SLAC accelerator department and the technical staffs of our collaborating institutions for their outstanding efforts on our behalf. We thank A. Brandenburg, P. Nason and C. Oleari for helpful discussions.

Work supported by Department of Energy contracts:

DE-FG02-91ER40676 (BU), DE-FG03-91ER40618 (UCSB), DE-FG03-92ER40689 (UCSC),

DE-FG03-93ER40788 (CSU), DE-FG02-91ER40672 (Colorado), DE-FG02-91ER40677 (Illinois),

DE-AC03-76SF00098 (LBL), DE-FG02-92ER40715 (Massachusetts), DE-FC02-94ER40818 (MIT),

DE-FG03-96ER40969 (Oregon), DE-AC03-76SF00515 (SLAC), DE-FG05-91ER40627 (Tennessee),

DE-FG02-95ER40896 (Wisconsin), DE-FG02-92ER40704 (Yale);

National Science Foundation grants:

PHY-91-13428 (UCSC), PHY-89-21320 (Columbia), PHY-92-04239 (Cincinnati),

PHY-95-10439 (Rutgers), PHY-88-19316 (Vanderbilt), PHY-92-03212 (Washington);

The UK Particle Physics and Astronomy Research Council (Bristol, Brunel, QMUL, Oxford, RAL);

The Istituto Nazionale di Fisica Nucleare of Italy

(Bologna, Ferrara, Frascati, Pisa, Padova, Perugia);

The Japan-US Cooperative Research Project on High Energy Physics (Nagoya, Tohoku);

The Korea Research Foundation (Soongsil, 1997).

## ∗∗List of Authors

Koya Abe,(25) Kenji Abe,(15) T. Abe,(22) I. Adam,(22) H. Akimoto,(22) D. Aston,(22) K.G. Baird,(11) C. Baltay,(31) H.R. Band,(30) T.L. Barklow,(22) J.M. Bauer,(12) G. Bellodi,(17) R. Berger,(22) G. Blaylock,(11) J.R. Bogart,(22) G.R. Bower,(22) J.E. Brau,(16) M. Breidenbach,(22) W.M. Bugg,(24) D. Burke,(22) T.H. Burnett,(29) P.N. Burrows,(19) A. Calcaterra,(8) R. Cassell,(22) A. Chou,(22) H.O. Cohn,(24) J.A. Coller,(4) M.R. Convery,(22) V. Cook,(29) R.F. Cowan,(13) G. Crawford,(22) C.J.S. Damerell,(20) M. Daoudi,(22) N. de Groot,(2) R. de Sangro,(8) D.N. Dong,(13) M. Doser,(22) R. Dubois, I. Erofeeva,(14) V. Eschenburg,(12) E. Etzion,(30) S. Fahey,(5) D. Falciai,(8) J.P. Fernandez,(27) K. Flood,(11) R. Frey,(16) E.L. Hart,(24) K. Hasuko,(25) S.S. Hertzbach,(11) M.E. Huffer,(22) X. Huynh,(22) M. Iwasaki,(16) D.J. Jackson,(20) P. Jacques,(21) J.A. Jaros,(22) Z.Y. Jiang,(22) A.S. Johnson,(22) J.R. Johnson,(30) R. Kajikawa,(15) M. Kalelkar,(21) H.J. Kang,(21) R.R. Kofler,(11) R.S. Kroeger,(12) M. Langston,(16) D.W.G. Leith,(22) V. Lia,(13) C. Lin,(11) G. Mancinelli,(21) S. Manly,(31) G. Mantovani,(18) T.W. Markiewicz,(22) T. Maruyama,(22) A.K. McKemey,(3) R. Messner,(22) K.C. Moffeit,(22) T.B. Moore,(31) M. Morii,(22) D. Muller,(22) V. Murzin,(14) S. Narita,(25) U. Nauenberg,(5) H. Neal,(31) G. Nesom,(17) N. Oishi,(15) D. Onoprienko,(24) L.S. Osborne,(13) R.S. Panvini,(28) C.H. Park,(23) I. Peruzzi,(8) M. Piccolo,(8) L. Piemontese,(7) R.J. Plano,(21) R. Prepost,(30) C.Y. Prescott,(22) B.N. Ratcliff,(22) J. Reidy,(12) P.L. Reinertsen,(27) L.S. Rochester,(22) P.C. Rowson,(22) J.J. Russell,(22) O.H. Saxton,(22) T. Schalk,(27) B.A. Schumm,(27) J. Schwiening,(22) V.V. Serbo,(22) G. Shapiro,(10) N.B. Sinev,(16) J.A. Snyder,(31) H. Staengle,(6) A. Stahl,(22) P. Stamer,(21) H. Steiner,(10) D. Su,(22) F. Suekane,(25) A. Sugiyama,(15) A. Suzuki,(15) M. Swartz,(9) F.E. Taylor,(13) J. Thom,(22) E. Torrence,(13) T. Usher,(22) J. Va’vra,(22) R. Verdier,(13) D.L. Wagner,(5) A.P. Waite,(22) S. Walston,(16) A.W. Weidemann,(24) E.R. Weiss,(29) J.S. Whitaker,(4) S.H. Williams,(22) S. Willocq,(11) R.J. Wilson,(6) W.J. Wisniewski,(22) J.L. Wittlin,(11) M. Woods,(22) T.R. Wright,(30) R.K. Yamamoto,(13) J. Yashima,(25) S.J. Yellin,(26) C.C. Young,(22) H. Yuta.(1)

(1)Aomori University, Aomori, 030 Japan, (2)University of Bristol, Bristol, United Kingdom, (3)Brunel University, Uxbridge, Middlesex, UB8 3PH United Kingdom, (4)Boston University, Boston, Massachusetts 02215, (5)University of Colorado, Boulder, Colorado 80309, (6)Colorado State University, Ft. Collins, Colorado 80523, (7)INFN Sezione di Ferrara and Universita di Ferrara, I-44100 Ferrara, Italy, (8)INFN Laboratori Nazionali di Frascati, I-00044 Frascati, Italy, (9)Johns Hopkins University, Baltimore, Maryland 21218-2686, (10)Lawrence Berkeley Laboratory, University of California, Berkeley, California 94720, (11)University of Massachusetts, Amherst, Massachusetts 01003, (12)University of Mississippi, University, Mississippi 38677, (13)Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, (14)Institute of Nuclear Physics, Moscow State University, 119899 Moscow, Russia, (15)Nagoya University, Chikusa-ku, Nagoya, 464 Japan, (16)University of Oregon, Eugene, Oregon 97403, (17)Oxford University, Oxford, OX1 3RH, United Kingdom, (18)INFN Sezione di Perugia and Universita di Perugia, I-06100 Perugia, Italy, (19)Queen Mary, University of London, London, E1 4NS, United Kingdom, (20)Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OX11 0QX United Kingdom, (21)Rutgers University, Piscataway, New Jersey 08855, (22)Stanford Linear Accelerator Center, Stanford University, Stanford, California 94309, (23)Soongsil University, Seoul, Korea 156-743, (24)University of Tennessee, Knoxville, Tennessee 37996, (25)Tohoku University, Sendai, 980 Japan, (26)University of California at Santa Barbara, Santa Barbara, California 93106, (27)University of California at Santa Cruz, Santa Cruz, California 95064, (28)Vanderbilt University, Nashville,Tennessee 37235, (29)University of Washington, Seattle, Washington 98105, (30)University of Wisconsin, Madison,Wisconsin 53706, (31)Yale University, New Haven, Connecticut 06511.


text:
# Properties of CO Molecular Gas in IR Luminous Galaxies

## 1. Introduction

Understanding the properties and evolution of the gas and dust content in nearby IR luminous galaxies is essential for understanding the cause and temporal evolution of starburst activity and its role in the cosmic evolution of galaxies. Therefore, there is a need to investigate large statistical samples of IR luminous galaxies using a multitude of different types of data, including CO, H I, and continuum in the submillimeter (sub-mm), FIR and radio, in order to constrain theories of how the interstellar medium (ISM) evolves. Millimeter and submillimeter CO rotational lines are often used as tracers of molecular hydrogen. The ratio of 12CO(3-2) to (1-0) line emission provides more sensitive measure of gas temperature and density than the ratio of 12CO(2-1) to (1-0) lines. Recently, we have presented the largest statistical CO survey (Yao et al. 2003) for the nearby universe by investigating the properties of CO molecular gas in a nearly complete subsample of 60 IR luminous galaxies selected from the SLUGS Survey (Dunne et al. 2000). The SLUGS survey containing 104 galaxies is based on the Revised Bright Galaxy Sample of IRAS galaxies (Soifer et al. 1987) within -10 $``$ $`\delta `$ $``$ +50 and with $`cz`$ $``$ 1900 km s-1, a flux limit of $`S_{60\mu m}`$ $``$ 5.24 Jy, and the FIR luminosity $`L_{FIR}`$ $``$ 1010 L.

## 2. Properties of Molecular Gas and Dust

A study of molecular gas of a complete flux-limited subsample is important to complement the study of the dust. The angular resolution of the 12CO(1-0) and (3-2) point observations presented in our CO survey are nearly identical ($``$ 15). A complete discussion of the observations, the CO line spectra, and the data analysis is contained in Yao et al. (2003) and references therein. Together with the sub-mm fluxes (Dunne et al. 2000), plus existing data on H I (Thomas et al. 2002), radio continuum (Condon et al. 1990), and FIR (Dunne et al. 2000; NED111The NASA/IPAC Extragalactic Database operated by the Jet propulsion Laboratory.), we are able to search for a relationship between the degree of excitation of the CO in this sample and the star forming properties.

### 2.1. CO and FIR Luminosities

Figure 1(a) shows a comparison between the CO and FIR luminosities for two different FIR color regimes. A linear fit to the 1-0 data in the log-log plane yields log$`{}_{10}{}^{}L_{{}_{}{}^{12}CO(10)}^{}`$ = (0.2 $`\pm `$ 0.7) + (0.57 $`\pm `$ 0.07) log$`{}_{10}{}^{}L_{FIR}^{}`$. The nonlinear variation of $`L_{CO}`$ with $`L_{FIR}`$ implies a decrease in the ratio $`L_{CO}`$/$`L_{FIR}`$ with increasing $`L_{FIR}`$, so that more IR luminous galaxies have higher dust temperatures and star formation efficiency. No correlation with the projected beam size on the galaxies is evident. The ratio depends intrinsically on the total FIR luminosity. This implies that for high-$`z`$ objects where $`L_{FIR}`$ is exceedingly high, the corresponding $`L_{CO}`$ will be seriously overestimated if a linear relation between $`L_{CO}`$ and $`L_{FIR}`$ is assumed.

### 2.2. Star Formation Efficiency

We examined the excitation of CO and its relation with the properties of gas/dust and star formation in the the central starburst regions in the SLUGS sample. There are no significant correlations between the 12CO(3-2) to (1-0) line ratio $`r_{31}`$ and the distance of galaxies, star formation rate, $`T_{dust}`$ and $`M_{dust}`$, H2 gas mass, and the color indices. This might reflect a range of localized conditions in the molecular clouds. But the degree of CO excitation measured by $`r_{31}`$ varies greatly from galaxy to galaxy. There is a trend for the $`r_{31}`$ to increase with increasing concentration and efficiency of star-forming activity (see Figure 1 (b)). The saturation of $`r_{31}`$ seen at higher SFE implies that the gas is denser and warmer in regions of higher SFE, which is consistent with the nonlinearity between $`L_{CO}`$ and $`L_{FIR}`$.

### 2.3. The CO-to-H2 Conversion Factor $`X`$

Our CO line measurements, together with dust and CO isotope data taken from the literature (Dunne et al. 2000; Aalto et al. 1995; Taniguchi et al. 1999), are modeled using the LVG approximation yielding $`X`$ = $`\frac{n(H_2)\mathrm{\Lambda }}{X_{CO}T_{rad}}`$ (see Figure 1 (c) and (d)) to estimate that the controversial $`X`$ lies between 1.3$`\times `$1019 and 6.7$`\times `$1019 cm-2(K km s-1)-1 for SLUGS galaxies, which is about 4-20 times lower than the conventional $`X`$ factor derived from the Galaxy. The mean value 2.7$`\times `$1019 is comparable to that estimated from diffuse clouds in the Galaxy and with that found for extreme starbursts in nearby galaxies.

### 2.4. Virial Stability of the Molecular Clouds

Using virial analysis, we show that the molecular clouds in starburst regions are not gravitationally bound, confirming the suggestion recently made by Zhu et al. (2003) for the Antennae galaxies, unless one is willing to accept a 9-90 times lower \[CO/H2\] abundance ratio. The cause could be the destruction by the stellar winds and expanding shells of newborn massive stars. Most of the CO line emission originates from the nonvirialized warm and diffuse gas clouds.

### 2.5. Revised Masses of Dust and Gas of SLUGS Galaxies

Our 12CO(3-2) observations provide an important database for correcting the SCUBA 850 $`\mu `$m data in the SLUGS sample (see Figure 1 (e) and (f)), thus permitting a revised characterization of the masses of dust and gas (see Seaquist et al. 2003). The average downward correction of dust mass is 25-38%, which has no bearing on earlier conclusions concerning the shapes of the dust mass-luminosity function derived from SLUGS survey. Using the 12CO(3-2)/(1-0) ratios in Yao et al. (2003), we estimate that it is unlikely that any correction is required for the contribution by 12CO(6-5) to 450 $`\mu `$m fluxes of the SLUGS galaxies measured with SCUBA. The revised mean gas-to-dust ratio (g/d) is reduced from 310 to 100, which is significantly lower than those reported for the global g/d ratios for IR luminous galaxies.

#### Acknowledgments.

I thank my advisor Prof. Ernie Seaquist and Reinhardt Travel Fellowship at Univ. of Toronto for their support. I also thank Kavli Institute for Theoretical Physics for their generosity.

## References

Aalto, S., et al. 1995, A&A, 300, 369

Condon, J. J., Helou, G., Sanders, D. B., & Soifer, B. T. 1990, ApJS, 73, 359

Dunne, L., et al. 2000, MNRAS, 315, 115

Seaquist, E. R., Yao, Lihong, Dunne, L., & Cameron, H. 2003, MNRAS, in press

Soifer, B. T., et al. 1987, ApJ, 320, 238

Taniguchi, Y., Ohyama, Y., & Sanders, D. B. 1999, ApJ, 522, 214

Thomas, H. C., et al. 2002, MNRAS, 329, 747

Yao, Lihong, Seaquist, E. R., Kuno, N., & Dunne, L. 2003, ApJ, 588, 771

Zhu, M., Seaquist, E. R., & Kuno, N. 2003, ApJ, 588, 243