text:
# Mixing patterns and community structure in networks

## 1 Introduction

Much of the recent research on the structure of networks of various kinds has looked at properties like path lengths, transitivity, degree distributions, and resilience of networks to vertex deletion , all of which, while of exceptional importance in many contexts, tend to focus our attention on the properties of individual vertices or vertex pairs—how far apart they are, what their degrees are, and so forth. However, in other contexts it may be equally important to ask about the large-scale properties of the network as a whole. Numbers of components and their distribution of sizes would be an example of such a property, one which is relevant to issues of accessibility  and to epidemiology . Searchability and the performance of search algorithms on networks would be another . A third is the existence and effects of large-scale inhomogeneity in networks—what we call “community structure”, the presence (or absence) in the network of regions with high densities of connections between vertices and other regions with low densities—and it is with a discussion of this topic that we begin this paper. (In some circles, this phenomenon is called “clustering”, an unfortunate terminology which risks confusion with another use of the word clustering introduced recently by Watts and Strogatz . We will use the word clustering only in reference to hierarchical clustering, which is a standard technique for community detection; otherwise we will avoid it.) Our investigation of community structure will lead us to consideration of mixing patterns in networks—which vertices connect to which others and why—as an explanation for observed communities in networks of all kinds, and eventually to consideration of more general classes of correlated networks including networks with correlations between the degrees of adjacent vertices.

Much of the work reported in this article has appeared previously in various papers, which the reader may like to consult for more detail than we can give here .

## 2 Community structure

The oldest studies by far of the large-scale statistical properties of networks are the studies of social networks carried out within the sociological community, which stretch back at least to the 1930s . Social networks are network representations of relationships of some kind, generically called “ties”, between people or groups of people, generically called “actors”. Actors might be individuals, organizations or companies, while ties might represent friendship, acquaintance, business relationships or financial transactions, amongst other things.

A long-standing goal among social network analysts has been to find ways of analysing network data to reveal the structure of the underlying communities that they represent. It is commonly supposed that the actors in most social networks divide themselves naturally into groups of some kind, such that the density of ties within groups is higher than the density of ties between them. A sketch of a network with such community structure is shown in Fig. 1.

It is a matter of common experience that social networks do contain communities. We look around ourselves and see that we belong to this clique or that, that we have a circle of close friends and others whom we know less well, that there are groupings within our personal networks on the basis of interest, occupation, geographical location and so forth. This does not however guarantee that a network contains community structure of type that we are considering here. It would be perfectly possible for each person in a network to have a well-defined set of close acquaintances, their own personal network neighbourhood, but for the network neighbourhoods of different people to overlap only partially, so that the network as a whole is quite homogeneous, with no clear communities emerging from the pattern of vertices and edges. A network model showing precisely this type of structure has been proposed and studied recently by Kleinberg . Our purpose in this section will be to investigate methods for detecting whether true community structure does exist in networks and for extracting the communities, and to apply those methods to particular networks. As we will see, the early intuition of the sociologists was correct, and many of the networks studied, including non-social networks, do possess large-scale inhomogeneity of precisely the type that would indicate the presence of community divisions.

The problem then is to take a network, specified in the simplest case by a list of $`n`$ vertices joined in pairs by $`m`$ edges, and from this structure to extract a set of communities—non-overlapping subsets of vertices that are, in some sense, tightly knit, having stronger within-group connections than between-group connections. The traditional, and still most common, method for detecting structure of this kind is the method of “hierarchical clustering” . In this method one defines a connection strength for each pair of vertices in the network, i.e., $`\frac{1}{2}n(n1)`$ numbers that represent a distance or weight for the connection between each pair. (In some versions of the method not all pairs are assigned a connection strength, in which case those that are not can be assumed to have a connection strength of zero.) Examples of possible definitions for the strengths include geodesic (shortest path) distances between pairs, or their inverses if one wants a measure that increases when pairs are more closely connected, counts of numbers of vertex- or edge-independent paths between pairs (“maxflow” methods) or weighted counts of total numbers of paths between pairs (adjacency matrix methods).

Then, starting with the $`n`$ vertices but no edges between them, one joins vertices together in decreasing order of the weights of vertex pairs, ignoring the edges of the original network. One can pause at any stage in this process and observe the pattern of components formed by the connections added so far, which are taken to be the communities of the network at that stage. The heirarchical clustering method thus defines not just a single decomposition of the network into communities, but a nested hierarchy of possible decompositions, having varying numbers of communities. This hierarchy can be represented as a tree or “dendrogram”, an example of which is shown in Fig. 2. A horizontal cut through the dendrogram at any given height, such as that denoted by the dotted line in Fig. 2, splits the tree into the communities for the corresponding stage in the hierarchical clustering process. By varying the height of the cut, one can arrange for the number of communities to take any desired value.

The construction of dendrograms is a popular technique for the analysis of network data, particularly within the sociological community. Software packages for network analysis, such as Pajek and UCInet, incorporate hierarchical clustering as a standard feature: for any network one can calculate a huge variety of vertex–vertex weights of different types and construct the corresponding dendrogram for any of them. The method however has some problems. There are many cases in which networks have rather obvious community structure, but hierarchical clustering fails to find it. One particular pathology that is frequently observed is that peripheral vertices tend to get disconnected from the bulk of the network, rather than being associated with the groups or communities that they are primarily attached to. For example, if a vertex is connected to the rest of the network by only a single edge, then presumably, were one to assign it to a community, it would be assigned to the community that the single edge leads to. In many cases, however, the hierarchical clustering method will declare the vertex instead to be a single-vertex community in its own right, in disagreement with our intuitive ideas of community structure.

In a recent paper therefore  we have proposed an alternative method for detecting community structure, based on calculations of so-called edge betweenness for vertex pairs. As we will see, this method detects the known community structure in a number of networks with remarkable accuracy.

### 2.1 Edge betweenness and community detection

Freeman  proposed a measure of centrality for the actors in a social network which he called “betweenness”. The betweenness of an actor is defined to be the number of shortest paths between pairs of vertices that pass through that actor. In cases where the number $`p`$ of shortest paths between a vertex pair is greater than one, each path is given an equal weight of $`1/p`$. Trivial algorithms for calculating betweenness take $`\mathrm{O}(mn^2)`$ time to calculate betweenness for all vertices, or $`\mathrm{O}(n^3)`$ time on a sparse graph (i.e., one in which the number of edges per vertex is constant in the limit of large graph size). This makes the calculation prohibitively costly on large networks. Recently however, two new algorithms have been proposed  that both allow the same calculation to be performed faster, in time $`\mathrm{O}(mn)`$, or $`\mathrm{O}(n^2)`$ on a sparse graph, by eliminating needless recalculations of geodesic paths. The betweenness of a vertex gives an indication, as the name implies, of how much the vertex is “between” other vertices. If, for example, information (or anything else) spreads through a network primarily by following shortest paths, then betweenness scores will indicate through which vertices most information will flow on average. The vertices with highest betweenness are also those whose removal will result in an increase to the geodesic distance between the largest number of other vertex pairs.

Here we consider an extension of Freeman’s betweenness to the edges in a network. The betweenness of an edge is defined to be the number of shortest paths between pairs of vertices that run along that edge, with paths again being given weights $`1/p`$ when there are $`p>1`$ between a given pair of vertices. In fact, the concept of edge betweenness actually appears to predate Freeman’s work on vertex betweenness, having appeared in an obscure technical report by an Amsterdam mathematician some years earlier . Edge betweenness has received very little attention in other literature until recently, but it provides us with an excellent measure of which edges in a network lie between different communities. In a network with strong community structure—groups of vertices with only a few inter-group edges joining them—at least some of the inter-group edges will necessarily receive high edge betweenness scores, since they must carry the geodesic paths between vertex pairs that lie in different communities. This implies that eliminating edges with high edge betweenness from a graph will remove the inter-group edges, and hence split the graph efficiently into its different groups. This is the principle behind our method for the detection of community structure. Our algorithm is as follows.

1. We calculate the edge betweenness of every edge in the network.
2. We remove the edge with the highest betweenness score, or randomly choose one such if more than one edge ties for the honour.
3. We recalculate betweenness scores on the resulting network and repeat from step 2 until no edges remain.

The recalculation in step 3 is crucial to the method’s success. When there is more than one inter-group edge between two groups of vertices, there is no guarantee that both will receive high betweenness scores; in some cases most geodesic paths will flow along one edge and only that one will receive a high score. Recalculation ensures that at some stage in the working of the algorithm each inter-group edge receives a high score and thus gets removed.

The calculation of all edge betweennesses takes time $`\mathrm{O}(mn)`$, and its repetition for all $`m`$ edges thus gives the algorithm a worst-case running time of $`\mathrm{O}(m^2n)`$, or $`\mathrm{O}(n^3)`$ on a sparse graph. The results of the algorithm can be represented as a dendrogram, just as in traditional hierarchical clustering, although one should be aware that the construction of the tree is not logically the same: the recalculation of the betweennesses after each edge removal means that there is no single function that can be defined for each edge in the initial graph such that the resulting dendrogram is the representation of a hierarchical clustering construction carried out using that function.

### 2.2 Examples

Here we give three examples of the application of our community structure finding algorithm to different networks. The first example is a set of computer generated graphs, specifically created to test the algorithm. We created a large number of graphs of 128 vertices each, divided into four groups of 32 vertices. Edges were placed at random between vertices within the same group with probability $`p_{\mathrm{in}}`$ and between vertices in different groups with probability $`p_{\mathrm{out}}`$, with the values of $`p_{\mathrm{in}}`$ and $`p_{\mathrm{out}}`$ chosen to make the average degree of a vertex equal to 16, and $`p_{\mathrm{out}}p_{\mathrm{in}}`$. These graphs were then fed into our community structure algorithm, and we measured what fraction of the vertices were correctly classified into their communities as a function of the ratio of $`p_{\mathrm{in}}`$ to $`p_{\mathrm{out}}`$, or equivalently the mean number $`z_{\mathrm{out}}`$ of edges from a vertex to vertices in other communities. The results are shown in Fig. 3. As the figure shows, the algorithm performs almost perfectly for values of $`z_{\mathrm{out}}`$ up to about 6. Beyond this point, as $`z_{\mathrm{out}}`$ approaches the value of 8 at which each vertex has as many inter-group edges as intra-group ones, the fraction of successfully classified vertices falls off sharply.

On the same plot we also show the performance of a standard hierarchical clustering algorithm based on edge-independent path counts (maxflow) on the same set of random graphs. As the figure shows, the traditional method is far inferior to our new algorithm in finding the known community structure.

For our second example, we move to real-world network data. In 1977, Wayne Zachary published the results of an ethnographic study he had conducted of social interactions between 34 members of a karate club at an American university . He recorded social contacts between members of the club over a two year period and published his results in the form of social networks. Fortuitously there arose, during the course of the study, a dispute between the two leaders of the club, the karate teacher and the club’s president, over whether to raise the club’s fees. Ultimately, the dispute resulted in the departure of the karate teacher and his starting another club of his own, taking with him about a half of the original club’s members. Here we analyse a network constructed by Zachary of friendships between club members before the split occurred. We compare the predictions of our community-finding algorithm applied to this network with the known lines along which the club divided. Our results are shown in Fig. 4.

In panel (a) of the figure we show the original network, with the grey squares representing the faction that ultimately sided with the teacher (who is vertex number 1), and the open circles the faction that sided with the club’s president (vertex number 34). In panel (b) we show the dendrogram output by our algorithm for this network. As the figure shows, the algorithm again performs nearly perfectly, with only one vertex, vertex number 3, being misclassified. (Inspection of panel (a) reveals that vertex 3 is in fact precisely caught in the middle of the network between the two factions, and so it is not entirely surprising that this vertex was misclassified.) Bear in mind that the network in this example was recorded *before* the fission of the club, so that the results of panel (b) are in some sense a prediction of events that were, at that time, yet to occur.

Finally, for our third example, we take a network for which we do not have any strong presuppositions about a “correct” division into communities. This example is a true experiment to see what information the algorithm can give us about a network whose structure is not wholly understood. The network in question is a food web, the web of trophic interactions (who eats whom) of marine organisms living in the Chesapeake Bay. The network was assembled by Baird and Ulanowicz  and contains 33 vertices representing the ecosystem’s most prominent taxa. The edges in a food web are, technically, directed; they can be thought of as pointing from prey to their predators, thus indicating the direction of energy (or carbon) flow up the food chain. Here however we have ignored the directed nature of the network, considering the edges merely to be undirected indicators of trophic interaction between taxon pairs.

The dendrogram produced for this food web by our community structure algorithm is shown in Fig. 5. As we can see, the algorithm splits the network into two principle communities and a couple of smaller peripheral ones. We have shaded the vertices in the dendrogram to show which taxa are surface dwellers in the bay (pelagic species) and which bottom dwellers (benthic species). A few species are of undetermined status. It is clear that our algorithm has in this case primarily extracted from the network the distinction between pelagic and benthic taxa. Thus our results appear to imply that the food web in question can be split roughly into separate surface- and bottom-dwelling subsystems, with relatively weak interaction between the two. A small number of benthic species are found to belong more strongly to the pelagic community than to the benthic one, perhaps indicating that a simple classification of species by where they live is not telling the whole story for this system. The results of our analysis might also be helpful in assigning a type to the undetermined species in the network.

## 3 Origins of community structure and assortative mixing

There is certainly more than one possible explanation for the presence of community structure in a network, and different explanations may be appropriate to different networks. In the case of a social network, for example, Jin et al.  have shown that communities can arise as a result of the growth dynamics of a network. If an acquaintance network grows by the introduction of pairs of people to one another by a mutual acquaintance, then an initial chance acquaintance with a member of a certain community will lead to introductions to other members of that community, so that one ultimately becomes linked to many of the community’s members and so becomes a member oneself. Using a simple computer model of this process, Jin et al. found that even networks with no initial community structure quickly develop such structure over time. One can think of this as a mechanism for the development of cliquishness in social networks.

This mechanism however is quite specific to social networks and could not be easily applied, for example, to the food web studied in the last section. It also completely ignores any personal attributes of the actors involved or affinities between actor pairs. A more general and perhaps more convincing explanation for community formation, which takes these things into account, is that of assortative mixing,111The name “assortative mixing” comes from the epidemiology community, where this effect has been studied extensively. It is also sometimes called “assortative matching”, particularly by ecologists. which is the tendency for nodes in a network to form connections preferentially to others that are like them in some way.

An example of assortative mixing in social networks is mixing by race. Table 1 shows data from the AIDS in Multiethnic Neighborhoods study , on mixing by race among sexual partners in the city of San Francisco, California. This part of the study focused on heterosexual partnerships, and the rows and columns of the matrix represent men and women in such partnerships, grouped by their (self-identified) race. Diagonal elements of the matrix represent the fraction of survey respondents in partnerships with members of their own group, and off-diagonal those in partnerships with members of other groups. Inspection of the figures shows that the matrix has considerably more weight along its diagonal than off it, indicating that assortative mixing does take place in this network. One might well expect mixing of this type to result in divisions within the community along lines of race, and we will show shortly that, within the context of simulations of network formation, assortative mixing can indeed give rise to community structure.

The amount of assortative mixing in a network can be characterized by measuring how much of the weight in the mixing matrix falls on the diagonal, and how much off it. Let us define $`e_{ij}`$ to be the fraction of all edges in a network that join a vertex of type $`i`$ to a vertex of type $`j`$. In the case of the matrix of Table 1, where the ends of an edge always attach to one man and one woman, we also specify which index corresponds to which type of end, which makes $`e_{ij}`$ asymmetric. For example, we could specify that the first index $`i`$ represents the man and the second $`j`$ the woman. For networks in which there is no corresponding distinction, $`e_{ij}`$ will be symmetric. The matrix should also satisfy the sum rules

$$\underset{ij}{}e_{ij}=1,\underset{j}{}e_{ij}=a_i,\underset{i}{}e_{ij}=b_j,$$
(1)
where $`a_i`$ and $`b_i`$ are the fraction of each type of end of an edge that is attached to vertices of type $`i`$. The values of $`a_i`$ and $`b_i`$ for the San Francisco study are also shown in Table 1. On graphs where there is no distinction between the ends of edges, we will have $`a_i=b_i`$.

Now we can define a quantitative measure $`r`$ of the level of assortative mixing in the network thus :

$$r=\frac{\underset{i}{}e_{ii}\underset{i}{}a_ib_i}{1_ia_ib_i}=\frac{Tr𝐞𝐞^2}{1𝐞^2},$$
(2)
where $`𝐞`$ is the matrix whose elements are the $`e_{ij}`$, and the notation $`𝐱`$ indicates the sum of the elements of the matrix $`𝐱`$. We call the quantity $`r`$ the “assortativity coefficient”. It takes the value 1 in a perfectly assortative network, since in that case the entire weight of the matrix $`𝐞`$ lies along its diagonal and $`_ie_{ii}=1`$. Conversely, if there is no assortative mixing at all, then $`e_{ij}=a_ib_j`$ for all $`i,j`$ and $`r=0`$. Networks can also be disassortative: vertices may associate preferentially with others of different types—the “opposites attract” phenomenon. In that case, $`r`$ will take a negative value.

One can certainly imagine that assortative mixing might apply in other types of networks as well. For example, we saw in Section 2.2 that a food web of marine organisms apparently divided into communities along lines of location—which species were surface dwellers (pelagic) and which bottom dwellers (benthic). It seems reasonable to hypothesize that the evolution of new predatory relationships between species is biased by the location of those species’ living quarters, and hence that the network structure would indeed reflect the pelagic/benthic division as a result of assortative mixing by location.

We can test our hypothesis that assortative mixing could be responsible for community formation in networks by computer simulation. Given a mixing matrix of the type shown in Table 1, we can create a random network with the corresponding mixing pattern and any desired degree distribution by the following algorithm.

1. First we choose degree distributions $`p_k^{(i)}`$ for each vertex type $`i`$. The quantity $`p_k^{(i)}`$ here denotes the probability that a randomly chosen vertex of type $`i`$ will have degree $`k`$. We can also calculate the mean degree $`z_i=_kkp_k^{(i)}`$ for each vertex type.
2. Next we choose a size for our graph in terms of the number $`m`$ of edges and draw $`m`$ edges from the desired distribution $`e_{ij}`$. We count the number of ends of edges of each type $`i`$, to give the sums $`m_i`$ of the degrees of vertices in each class, and we calculate the expected number $`n_i`$ of vertices of each type from $`n_i=m_i/z_i`$ (rounded to the nearest integer).
3. We draw $`n_i`$ vertices from the desired degree distribution $`p_k^{(i)}`$ for type $`i`$. Normally the degrees of these vertices will not sum exactly to $`m_i`$ as we want them to, in which case we choose one vertex at random, discard it, and draw another from the distribution $`p_k^{(i)}`$, repeating until the sum does equal $`m_i`$.
4. We pair up the $`m_i`$ ends of edges of type $`i`$ at random with the vertices we have generated, so that each vertex has the number of attached edges corresponding to its chosen degree.
5. We repeat from step 3 for each vertex type.

We have used this algorithm to generate example networks with desired levels of assortative mixing. For example, Fig. 6 shows an undirected network of $`n=100`$ vertices of four different types, generated using the symmetric mixing matrix

$$𝐞=\left(\begin{array}{cccc}0.18& 0.02& 0.01& 0.03\\ 0.02& 0.20& 0.03& 0.02\\ 0.01& 0.03& 0.16& 0.01\\ 0.03& 0.02& 0.01& 0.22\end{array}\right),$$
(3)
which gives a value of $`r=0.68`$ for the assortativity coefficient. A simple Poisson degree distribution with mean $`z=5`$ was used for all vertex types. The graph was then fed into the community finding algorithm of Section 2.1 and a cut through the resulting dendrogram performed at the four-community level. The communities found are shown by the four shapes of vertices in the figure and correspond very closely to the real vertex type designations, which are represented by the four different vertex shades. In other words, by introducing assortative mixing by vertex type into this network, we have created vertex-type communities that register in our community finding algorithm in exactly the same way as communities in naturally occurring networks. This strongly suggests that assortative mixing could indeed be an explanation for the occurrence of such communities, although it is worth repeating that other explanations are also possible.

## 4 Other types of assortative mixing

Assortative mixing can depend on vertex properties other than the simple enumerative properties discussed in the preceding section. For example, we can also have assortative mixing by scalar characteristics, either discrete or continuous. A classic example of such mixing, much studied in the sociological literature, is acquaintance matching by age. In many contexts, people appear to prefer to associate with others of approximately the same age as themselves. For example, consider Fig. 7, which shows the ages at marriage of the male and female members of 1141 married couples drawn from the US National Survey of Family Growth . Each point in the figure represents one couple, its position along the horizontal and vertical axes corresponding to the ages of the husband and wife respectively. The study was based on interviews with women, and was limited to those of childbearing age, so the vertical axis cuts off around 40. Also only the first marriage for each woman interviewed is shown, even if she married more than once. Despite these biases however, the figure reveals a clear trend: people prefer to marry others of an age close to their own.

It is perhaps stretching a point a little to consider first marriage ties between couples as forming a social network, since people have at most one first marriage and hence would have a maximum degree of one within the network. Here, however, we consider marriage age as a proxy for the ages of sexual partners in general, and conjecture that a similar age preference will be seen in non-married partners also. A recent study by Garfinkel et al.  appears to support this conjecture.

Assortative mixing according to scalar characteristics can result in the formation of communities, just as in the case of discrete characteristics. One could have separate communities formed of old and young people, for instance. However, it is also possible that we do not get well-defined communities, but instead get an overlapping set of groups with no clear boundaries, ranging for example from low age to high age. In the sociological literature such a continuous gradation of one community into another is called “stratification” of the network.

As with assortative mixing on discrete characteristics, one can define an assortativity coefficient to quantify the extent to which mixing is biased according to scalar vertex properties. To do this, we define $`e_{xy}`$ to be the fraction of edges in our network that connect a vertex of property $`x`$ to another of property $`y`$. The matrix $`e_{xy}`$ must satisfy sum rules as before, of the form

$$\underset{xy}{}e_{xy}=1,\underset{y}{}e_{xy}=a_x,\underset{x}{}e_{xy}=b_y,$$
(4)
where $`a_x`$ and $`b_y`$ are, respectively, the fraction of edges that start and end at vertices with ages $`x`$ and $`y`$. Then the appropriate definition for the assortativity coefficient is

$$r=\frac{\underset{xy}{}xy(e_{xy}a_xb_y)}{\sigma _a\sigma _b},$$
(5)
where $`\sigma _a`$ and $`\sigma _b`$ are the standard deviations of the distributions $`a_x`$ and $`b_y`$. The reader will no doubt recognize this definition of $`r`$ as the standard Pearson correlation coefficient for the quantities $`x`$ and $`y`$. It takes values in the range $`1r1`$ with $`r=1`$ indicating perfect assortative mixing, $`r=0`$ indicating no correlation between $`x`$ and $`y`$, and $`r=1`$ indicating perfect disassortative mixing, i.e., perfect anticorrelation between $`x`$ and $`y`$.

If we take the marriage data from Fig. 7, for example, and feed it into Eq. (5), we find that $`r=0.57`$, indicating once again that mixing is strongly assortative (as is in any case obvious from the figure).

Mixing could also depend on vector or even tensor characteristics of vertices. One example would be mixing by geographical location, which could be regarded as a two-vector. It seems highly likely that if one were to record both acquaintance patterns and geographical location for actors in a social network, one would discover that acquaintance is strongly dependent on geography, with people being more likely to know others who live in the same part of the world as themselves.

### 4.1 Mixing by vertex degree

We will spend the rest of this article examining one particular case of mixing according to a scalar vertex property, that of mixing by vertex degree, which has been studied for some time in the social networks literature and has recently attracted attention in the mathematical and physical literature also. Krapivsky and Redner  for instance found in studies of the preferential attachment model of Barabási and Albert  that edges did not fall between vertices independent of their degrees. Instead there was a higher probability to find some degree combinations at the ends of edges than others. Pastor-Satorras et al.  showed for data on the structure of the Internet at the level of autonomous systems that the degrees of adjacent vertices were anticorrelated, i.e., that high-degree vertices prefer to attach to low-degree vertices, rather than other high-degree ones—the network is disassortative by degree. To demonstrate this, they measured the mean degree degree $`k_{\mathrm{nn}}`$ of the nearest-neighbours of a vertex, as a function of that vertex’s degree $`k`$. They found that $`k_{\mathrm{nn}}`$ decreases with increasing $`k`$, approximately as $`k^{1/2}`$. That is, the mean degree of your neighbours goes down as yours goes up. Maslov and Sneppen  have offered an explanation of this result in terms of ensembles of graphs in which double edges between vertices are forbidden. Maslov and Sneppen also showed in a separate paper  that the protein interaction network of the yeast S. Cerevisiae displays a similar sort of disassortative mixing.

An alternative way to quantify assortative mixing by degree in a network is to use an assortativity coefficient of the type described in the previous section . Let us define $`e_{jk}`$ to be the fraction of edges in a network that connect a vertex of degree $`j`$ to a vertex of degree $`k`$. (As before, if the ends of an edge connect different types of vertices, then the matrix will be asymmetric, otherwise it will be symmetric.) In fact, we define $`j`$ and $`k`$ to be the “excess degrees” of the two vertices, i.e., the number of edges incident on them less the one edge that we are looking at at present. In other words, $`j`$ and $`k`$ are one less than the total degrees of the two vertices. This designation turns out to be mathematically convenient for many developments. If the degree distribution of the network as a whole is $`p_k`$, then the distribution of the excess degree of the vertex at the end of a randomly chosen edge is

$$q_k=\frac{(k+1)p_{k+1}}{z},$$
(6)
where $`z=_kkp_k`$ is the mean degree . Then one can define the assortativity coefficient to be

$$r=\frac{\underset{jk}{}jk(e_{jk}q_jq_k)}{\sigma _q^2},$$
(7)
where $`\sigma _q`$ is the standard deviation of the distribution $`q_k`$. On a directed or similar network, where the ends of an edge are not the same and $`e_{jk}`$ is asymmetric, this generalizes to

$$r=\frac{\underset{jk}{}jk(e_{jk}q_j^aq_k^b)}{\sigma _a\sigma _b},$$
(8)
where $`\sigma _a`$ and $`\sigma _b`$ are the standard deviations of the distributions $`q_k^a`$ and $`q_k^b`$ for the two types of ends. (The measure introduced by Pastor-Satorras et al.  can also be expressed simply in terms of the matrix $`e_{jk}`$: it is $`k_{\mathrm{nn}}=_jje_{jk}/q_k`$. Maslov and Sneppen  gave entire plots of the raw $`e_{jk}`$, using colours to code for different values. These plots are however rather difficult to interpret by eye.)

In Table 2 we show values of $`r`$ measured for a variety of different real-world networks. The networks shown are divided into social, technological and biological networks, and a particularly striking feature of the table is that the values of $`r`$ for the social networks are all positive, indicating assortative mixing by degree, while those for the technological and biological networks are all negative, indicating disassortative mixing. It is not clear at present why this should be, although explanations for the observed mixing behaviours have been proposed in some specific cases .

As with the mixing by discrete enumerative characteristics discussed in Section 3, we can also investigate the effects of assortative mixing by looking at computer generated networks with particular types of mixing. Unfortunately, no simple algorithm exists for generating graphs mixed by vertex degree analogous to that of Section 3 (see Refs.  and ) and one is forced to resort to Monte Carlo generation of graphs using Metropolis–Hastings type algorithms of the sort widely used for graph generation in mathematics and quantitative sociology. Such algorithms however are straightforward to implement. For the present case, we take the simple example form

$$e_{jk}=𝒩\mathrm{e}^{(j+k)/\kappa }\left[\left(\genfrac{}{}{0pt}{}{j+k}{j}\right)p^jq^k+\left(\genfrac{}{}{0pt}{}{j+k}{k}\right)p^kq^j\right],$$
(9)
where $`p+q=1`$, $`\kappa >0`$, and $`𝒩=\frac{1}{2}(1\mathrm{e}^{1/\kappa })`$ is a normalizing constant. This means that the distribution of the sum $`j+k`$ of the excess degrees at the ends of an edge falls off as a simple exponential, while that sum is distributed between the two ends binomially, the parameter $`p`$ controlling the assortative mixing. For values of $`p`$ ranging from 0 to $`\frac{1}{2}`$ we get various values of the assortativity $`r`$, both positive and negative, passing through zero at $`p_0=\frac{1}{2}\frac{1}{4}\sqrt{2}=0.1464\mathrm{}`$

As an example, we show in Fig. 8 the giant components of two graphs of this type generated using the Monte Carlo method. One of them, graph (a), is assortatively mixed by degree, while the other, graph (b), is disassortatively mixed. The difference between the two is clear to the eye. In the first case, because the high degree vertices prefer to attach to one another, there is a central “core” to the network, composed of these high-degree vertices, and a straggling periphery of low-degree vertices around it. In epidemiology a dense central portion of this type is called a “core group” and is thought to be capable of acting as a reservoir for disease, keeping diseases circulating even when the density of the network as a whole is too low to maintain endemic infection. In social network analysis one also talks of “core/periphery” distinctions in networks, another concept that mirrors what we see here. In the second graph, which is disassortative, a contrasting picture is evident: the high-degree vertices prefer not to associate with one another, and are as a result scattered widely over the network, producing a more uniform appearance.

To shed more light on the effects of assortativity, we show in Fig. 9 the size of the largest component in networks of this type as the degree distribution parameter $`\kappa `$ is varied, for various values of $`p`$. For low values of $`\kappa `$ the mean degree of the network is small, and the resulting density of edges is too low to produce percolation in the network, so there is no giant component. As $`\kappa `$ increases, however, there comes a point, clearly visible on the plot, at which the edge density is great enough to form a giant component. Figure 9 reveals two interesting features of this transition. First, the position of the transition, the value of the parameter $`\kappa `$ at which it takes place, is smaller in assortatively mixed networks than in disassortative ones. In other words, it appears that the presence of assortativity in the degree correlation pattern allows the network to percolate more easily. This result is intuitively reasonable: the core group of the assortative network seen in Fig. 8a has a higher density of edges than the network as a whole and so one would expect percolation to take place in this region before it would in a network with the same average density but no core group.

Second, the figure shows that, even though the assortative network percolates more easily than its disassortative counterpart, its largest component does not grow as large as that of the disassortative network in the limit where $`\kappa `$ becomes large. This too can be understood in simple terms: percolation occurs more easily when there is a core group, but is also largely confined to that core group and so does not spread to as large a portion of the network as it would in other cases.

In epidemiological terms, one could think of these two results as indicating that assortative networks will support the spread and persistence of a disease more easily than disassortative ones, because they possess a core group of connected high-degree vertices. But the disease is also restricted mostly to that core group. In a disassortative network, although percolation and hence epidemic disease requires a denser network to begin with, when it does happen it will affect a larger fraction of the network, because it is not restricted to a core group.

## 5 Conclusions

In this article we have examined two related properties of networks: community structure and assortative mixing. We have described a new algorithm for finding groups of tightly-knit vertices within networks—communities in our nomenclature—which is based on the calculation of an “edge betweenness” index for network edges. The algorithm appears to be successful at detecting known community structure in various example networks, and we have found that a number of real-world networks do indeed possess community structure to a greater or lesser degree.

Turning to possible explanations for this structure we have suggested that assortative mixing, the preferential association of vertices in a network with others that are like them in some way, is one possible mechanism for community formation. We have defined a measure of the strength of assortative mixing and applied it, for example, to data on mixing by race in social networks, showing that there is strong assortativity in this case, at least for the survey data that we have examined. We have also given a simple algorithm for creating networks with assortative mixing according to discrete characteristics imposed upon the vertices, and used it to generate example networks which, when fed into our community detection algorithm, reveal strong community structure similar to that seen in the real-world data. This lends some conviction to the theory that assortative mixing could, at least in some cases, be a contributing factor in the formation of communities within networks.

We have also looked at assortative mixing by scalar characteristics of vertices, such as the age of individuals in a social network, and particularly vertex degree. By measuring mixing of the latter type for a variety of different networks, we have shown that social networks appear often to be assortatively mixed by degree, while technological and biological networks appear normally to be disassortative. Using computer generated model networks we have also shown that assortativity by vertex degree makes networks percolate more easily—they develop a giant component for a lower average edge density than a similar network with neutral or disassortative mixing. Conversely, however, disassortative networks tend to have larger giant components when they do develop. These findings have implications for epidemiology, for example: they imply that a disease spreading on a network that is assortatively mixed, as most social networks appear to be, would reach epidemic proportions more easily than on a disassortative network, but that the epidemic might ultimately affect fewer people than in the disassortative case.

Looking ahead, some obvious next steps in the studies presented here are the application of community finding algorithms to other networks, the study of mixing patterns in other networks, and theoretical investigations of the effects of assortative mixing and other network correlations on network structure and function, including for instance network resilience and network epidemiology. A number of authors have already started work on these problems .

## Acknowledgements

The authors would like to thank Jennifer Dunne, Neo Martinez and Doug White for help assembling and interpreting the data used in Figs. 4 and 5, and László Barabási, Jerry Davis, Jennifer Dunne, Jerry Grossman, Hawoong Jeong, Neo Martinez and Duncan Watts for providing data used in the calculations for Table 2. The marriage data for Fig. 7 were provided by the Inter-University Consortium for Political and Social Research at the University of Michigan. This work was supported in part by the National Science Foundation under grants DMS–0109086 and DMS–0234188.


format: text

text:
# Tight binding studies of strained Ge/Si(001) growth

## 1 Introduction

Understanding the relationship between strain and the formation of nanostructures is extremely important as practical nanoelectronic devices are sought. The growth of Ge on Si(001) is particularly interesting in this context, not only because of the small, three dimensional features which form at high coverage (the so-called “hut clusters”) but also because of the direct compatibility with standard group IV semiconductor technology. In this paper, we present tight binding calculations of various aspects of Ge growth on Si(001) for coverages between one and three monolayers (ML) of Ge. This work follows naturally from our previous paper, where we studied the $`(2\times N)`$ reconstruction which forms at sub-monolayer coverages using first principles electronic structure techniques.

As the coverage of Ge on Si(001) approaches 1ML, it forms the well-known $`(2\times N)`$ reconstruction, where regularly spaced missing dimers appear in the surface; experimental observations of N show that it typically lies between 8 and 12, with the spacing dependent on growth source (gas source, e.g. GeH4, or solid source) and conditions (growth rate and temperature). As the thickness of the deposited layer increases, various further effects are seen on this reconstruction, again depending on growth source and conditions for their onset: the value of N decreases; the missing dimers align in neighbouring dimer rows to form straight trenches; a new, patch-like reconstruction with a double periodicity (called the $`(M\times N)`$ reconstruction) forms (illustrated in Fig. 1); the roughness of step edges changes (with the different step types eventually becoming equally rough, leading to steps running along the elastically soft (100) and (010) directions); and finally three dimensional structures such as “hut” pits and clusters form. We shall investigate the first three of these phenomena in this paper.

The modelling of the formation of straight trenches and the $`(M\times N)`$ reconstruction requires large unit cells, while the comparisons we are performing require a quantum mechanical technique for accuracy (for instance, to correctly model the dimer buckling on the surface). We have chosen to use an $`𝒪(N)`$ (or linear scaling) tight binding technique to meet these criteria (this is a technique where the computer effort required scales linearly with the system size at the expense of an approximation, rather than with the cube of the system size, as is the case for standard techniques, referred to as “diagonalisation”). However, this introduces two levels of approximation: the first, inherent in the tight binding parameterisation; the second, in the use of an $`𝒪(N)`$ method. We test these approximations separately. First, we test the parameterisation by comparing the results of exact diagonalisation to ab initio results. Second, we check the accuracy of the $`𝒪(N)`$ technique against the exact diagonalisation tight binding results. This gives us a clear idea of what errors are introduced and where. Starting from the $`(2\times N)`$ reconstruction, we investigate first the dependence of the value of N on the depth of the Ge layers. We then calculate the energetics of kinks in trenches of missing dimers with increasing Ge depth, in order to understand the formation of straight trenches of missing dimers. Finally we address the $`(M\times N)`$ reconstruction, calculating minimum energy values of M for a fixed value of N. In all of the above calculations, we will also address the question of an appropriate energy to use for Ge dimers when comparing reconstructions of different stoichiometry, which is closely related to the problem of finding an appropriate chemical potential for Ge.

The previous work in this area only addressed the $`(2\times N)`$ reconstruction, and investigated the dependence of N on depth of Ge using empirical potentials, such as the Stillinger-Weber and modified Keating forms . These studies found results in broad agreement with experiment, and our previous work provides support for the mechanisms suggested by these empirical studies. We note that intermixing of Si and Ge may well occur during experimental growth, though this depends on experimental conditions, and is thought to be largely suppressed after one monolayer of Ge. The likely effect would be to reduce the surface strain somewhat; we have chosen not to account for intermixing for three reasons. First, there are experimental results available for situations where intermixing is almost entirely suppressed, with which we can make direct contact. Second, it adds significant additional complication to the problem, where we are trying to find the underlying mechanisms. Third, the effects are likely to be relatively small. Indeed, in previous empirical work, the authors concluded that the Ge/Si interlayer mixing could not be substantial. In a recent detailed ab initio study of intermixing, it was shown that, while Ge can penetrate to the fourth layer, this is likely to be at the level of a few percent at a typical low growth temperature (e.g. 500C). Hence, we neglect this effect in the present work.

The rest of the paper proceeds as follows: in the next section we review the details of the computational techniques used, consider the question of what energies to use when comparing reconstructions of different periodicity (and hence stoichiometry) and test the accuracy of our methods; after this, we present our results for the different situations considered (the $`(2\times N)`$ reconstruction with increasing Ge depth; trench kinking with increasing Ge depth; the $`(M\times N)`$ reconstruction), and we finish with a discussion of the results.

## 2 Computational Details

### 2.1 Tight Binding Parameters

The tight binding technique has been discussed elsewhere in detail, as have the ideas and implementation of $`𝒪(N)`$ techniques. We used a parameterisation for Si-Si bonds (and Si-H bonds for termination) specifically designed for and tested extensively on the Si(001) surface, and a parameterisation for Ge-Ge and Ge-Si bonds again tested on the Si(001) surface, and fitted to reproduce the experimental Ge lattice constant (giving a lattice mismatch of 4% between Ge and Si). We used the Oxon code to perform exact diagonalisations, and the DensEl code (an implementation of the auxiliary density matrix (ADM) technique) for linear scaling calculations.

We work with periodic boundary conditions, simulating a surface via a slab with a vacuum gap equivalent to eight atomic layers of Si. In all cases, the slab modelled was ten atomic layers deep, with the bottom layer terminated in hydrogen. The length and width used depended on the calculations: for the $`(2\times N)`$ calculations, cells were one dimer row wide, and ranged from four to twelve dimers long (from 94 to 286 atoms); for trench kinking calculations, they were between two and six dimer rows wide and eight (1 or 2ML Ge) or six (3ML Ge) dimers long (from 284 atoms to 1,140 atoms); and for the $`(M\times N)`$ reconstructions, they were typically six dimer rows wide and eight dimers long (around 1,100 atoms). We investigated the importance of substrate relaxation, and found that in most cases (in all cases with 2 or more monolayers of Ge on Si(001)) allowing four layers of silicon below the Ge to relax gave energy convergence to within a few meV (and presented the same amount of substrate relaxation below the Ge). For the detailed comparison of one monolayer of Ge on Si(001) presented below, we found that allowing seven layers of silicon to relax was important. The reason for this lies in the rebonding of the second layer atoms across the gap where the dimer is missing. When these atoms are silicon, they are stretched across the gap, and require significant depth for relaxation of strain. When they are germanium, the rebonding helps to alleviate the mismatch with the silicon substrate, and less depth is required.

When using exact diagonalisation to test the parameterisation involving Ge against earlier DFT results, we used a $`4\times 2\times 1`$ Monkhorst-Pack k-point mesh, which we found converged total energies to within 2 or 3 meV. When performing the linear scaling calculations, we had to choose a cutoff on the density matrix, which in the DensEl code is typically specified as a number of hamiltonian-range hops. Previous work has shown that, for semiconductors and insulators, forces and energy differences are well converged at 3 hops (roughly 6 Å in Si) and completely converged at 5 hops (roughly 10 Å in Si). We present detailed calculations for the energy of one monolayer of Ge in the $`(2\times N)`$ reconstruction below (in section 2.3), but the conclusion is that 3 hops is an appropriate compromise between accuracy and CPU time.

### 2.2 Comparing different stoichiometries

In our previous paper, where we presented ab initio calculations of the $`(2\times N)`$ reconstruction for sub-monolayer coverages of Ge, we gave detailed arguments for the energy to be used for Ge dimers when comparing reconstructions of different periodicity and hence stoichiometry. This energy (which is closely related to a chemical potential) is extremely important when making contact with the experimentally observed values of N and M, which is one of the aims of this paper. We will briefly recall the statistical mechanical arguments used before, and then extend them to the rather different situations found in this work.

We are assuming that the surface is in thermal equilibrium, which means of course that we are ignoring kinetic effects — an assumption that will be discussed in Sec. 4. Then, for a given number of Ge atoms on the surface, the probability of finding them in a particular arrangement $`\mathrm{\Gamma }`$ will be proportional to $`\mathrm{exp}(E_\mathrm{\Gamma }/k_\mathrm{B}T)`$, with $`E_\mathrm{\Gamma }`$ the energy for the arrangement $`\mathrm{\Gamma }`$.

We only need to know how $`E_\mathrm{\Gamma }`$ varies from one arrangement to another since we are considering a fixed number of atoms. However, it will be helpful to use one specific arrangement as a “reference” arrangement; we choose this to be all Ge atoms paired into dimers and arranged in a perfect, lattice-matched monolayer covering a certain area of the surface. The shape is unimportant, but it is easier to consider a rectangle. We now consider forming regular arrays of missing dimer trenches by removing dimers from the monolayer and replacing them at the edges of the island, and allowing the system to relax. It is convenient to split this process into two parts: (i) fetch an appropriate number of Ge dimers (equivalent to the number to be removed) from infinity and place them at the boundary of the layer; (ii) remove the appropriate Ge dimers to form missing dimer trenches and relax the system. We now calculate the quantity $`\zeta (N)`$, which is the energy change per Ge dimer in going from a perfect layer to a periodically reconstructed one. It is given by:

$$\zeta (N)=\frac{(E_\mathrm{f}(N)+E_\mathrm{p})}{(N1)},$$
(1)
where $`E_\mathrm{f}(N)`$ is the formation energy of fully relaxed missing dimers in the top layer with spacing N (defined as the energy difference between a relaxed surface with a missing dimer reconstruction and a surface with a perfect top Ge layer), and $`E_\mathrm{p}`$ is the energy for a Ge dimer in the reference system (a perfect layer — the energy is defined as the difference per Ge dimer between a surface with a perfect top Ge layer, and the surface formed by removing the top layer of Ge and reconstructing and relaxing). By finding the minimum value of $`\zeta (N)`$ while varying N, we can find the expected value of N. It should be noted that if $`\zeta (N)`$ is positive, then the reconstruction is less stable than a perfect monolayer of Ge (which was chosen as the reference system).

Now we consider how to generalise this argument to the systems considered in this paper. There are four comparisons that we will make:

1. The $`(2\times N)`$ reconstruction, on multiple layers of Ge (compared to experiment) (Sec. 3.1)
2. The energy of missing dimer trenches with and without disorder (compared to each other) (Sec. 3.2)
3. The $`(M\times N)`$ reconstruction, with multiple layers of Ge (compared to experiment) (Sec. 3.3)
4. The energies of $`(2\times N)`$ and $`(M\times N)`$ reconstructions on the same number of layers of Ge (compared to each other) (Sec. 3.4)

For the second of these (the energetics of disorder in missing dimer trenches), we do not need to concern ourselves with compensation for different amounts of Ge, as there is no variation in period, just in arrangements of missing dimers. For the other three, we need to identify appropriate energies to use when making the comparisons discussed above, as well as the physical situation we are dealing with.

#### The $`(2\times N)`$ reconstruction on multiple layers of Ge

This is possibly the simplest case to deal with. To be clear, we have a system composed of a Si substrate, with $`\theta `$ layers of Ge on top, only the top layer of which has a missing dimer reconstruction (in other words we have $`\theta 1`$ full layers of Ge on top of Si as our substrate). Then we use the same argument as we used for sub-monolayer coverages of Ge on Si(001) — i.e. we consider a pool of Ge on top of $`\theta 1`$ full layers of Ge on Si, and examine different arrangements that we could make with it, calculating $`\zeta (N)`$ for each. The energy that we require for comparisons is then the energy difference per Ge dimer between $`\theta `$ full layers of Ge on Si, and $`\theta 1`$ full layers of Ge on Si, which we call $`E_\mathrm{p}(\theta )`$, where $`\theta `$ is the coverage (or number of layers).

#### The $`(M\times N)`$ reconstruction on multiple layers of Ge

The easiest energy to use for this case is little different to that used above; we assume that we have $`\theta 2`$ full layers of Ge on top of Si, and then a further layer of Ge with a $`(2\times N)`$ reconstruction (the $`\theta 1^{\mathrm{th}}`$ layer). This acts as our substrate; since there is a $`(2\times N)`$ reconstruction, and the dimer row direction rotates by $`90^{}`$ with each layer, we must form finite width strips of Ge on this substrate (where the width is $`N/21`$ dimer rows) to make an $`(M\times N)`$ reconstruction. We now consider our pool of Ge as divided into strips, and make different arrangements with that, now calculating $`\zeta (M)`$ for each. The energy required is then just the energy difference per Ge dimer between perfect strips of Ge on a $`(2\times N)`$ reconstructed surface with $`\theta 1`$ layers of Ge and the $`(2\times N)`$ surface reconstructed surface with $`\theta 1`$ layers of Ge, which we call $`E_{\mathrm{p}_{\mathrm{strip}}}(\theta )`$, where again $`\theta `$ is the coverage (or number of layers).

#### The $`(2\times N)`$ and $`(M\times N)`$ reconstructions on multiple layers of Ge

It is easiest to compare these reconstructions by comparing their stabilities relative to the perfect surface, which requires compensation for the missing Ge. For the $`(2\times N)`$ reconstruction, this requires only $`E_\mathrm{p}(\theta )`$, which has already been calculated above. For the $`(M\times N)`$ reconstruction, we need to perform a little more work, but only using quantities that we already have — the energy of Ge dimers in a perfect layer of a given thickness, $`E_\mathrm{p}(\theta )`$ (both $`\theta `$ and $`\theta 1`$ will be required: $`\theta 1`$ to compensate for the missing dimer trench in the underlying $`(2\times N)`$ reconstruction, and $`\theta `$ to compensate for the missing rows of dimers in the top layer), as well as $`E_{\mathrm{p}_{\mathrm{strip}}}(\theta )`$.

### 2.3 Testing Approximations

We now consider the errors introduced by the two approximations that we use: tight binding, rather than full ab initio theory; and the cutoff on the density matrix used in the linear scaling technique. In figure 2, we present the results for $`\zeta (N)`$ (defined above in equation 1) for sub-monolayer coverages of Ge on Si(001). The main graph shows the results found using exact diagonalisation plotted with the DFT results from our earlier paper, while the inset graph shows the $`𝒪(N)`$ results plotted with the exact diagonalisation results. The technique used to calculate $`\zeta (N)`$ is discussed in detail in Ref.  and briefly above in Sec. 2.2.

Considering the ab initio and tight binding exact diagonalisation results first, we can see that the agreement is good for large N, but that at $`N=4`$ in particular there is a considerable discrepancy. However, the reason for this is not due to the approximation made in using tight binding. The lattice mismatch of the ab initio results is too large (calculated to be 5.7% as opposed to the experimental and (fitted) tight binding value of 4%) which will increase any strain-related effects; also, the ab initio calculations allowed only the top four layers of the slab to relax. Using the same conditions as the ab initio calculations (top four layers relaxing, and rescaling the tight binding parameters to create a 5.7% lattice mismatch) we find $`\zeta (N=4)=0.054`$ eV, where the ab initio value is 0.079 eV, which is in much better agreement. This suggests that the discrepancy between the tight binding and ab initio results for small N is in large part due to the different lattice parameter and simulation conditions, and that the tight binding parameterisation has introduced only minor inaccuracies relative to the ab initio calculations.

We must also consider the effect of cutting off the density matrix (the essential approximation in creating an $`𝒪(N)`$ method) on the tight binding results. Data from the DensEl calculations for two different density matrix cutoffs (three and five hamiltonian range hops, which correspond roughly to 6 Å and 10 Å) and the exact diagonalisation are shown in the inset to Fig. 2. We can see that, at 5 hops, the $`𝒪(N)`$ results and the exact diagonalisation results are in almost exact agreement, while at 3 hops the overall shape of the curve is reproduced well, though the exact numerical values are a little high. However, these results are sufficiently accurate for our purposes, and we will use a 3 hop cutoff in all calculations in the paper. These results are in excellent agreement with previous work which suggested that energy differences were completely converged at 5 hops and well converged at 3 hops.

## 3 Results

In this section we present results for a series of increasingly complex reconstructions for coverages of Ge on Si(001) that exceed 1ML: calculations of the value of N for the $`(2\times N)`$ reconstruction; calculations of the energy of kinks in the missing-dimer trenches for this reconstruction (for 1ML as well as more than 1ML); and finally the $`(M\times N)`$ “patch” reconstruction.

### 3.1 The $`(2\times N)`$ reconstruction for more than 1ML

Under certain growth conditions, the $`(2\times N)`$ reconstruction is seen to persist when more than one monolayer of Ge has been deposited. This results in a more strained layer, and the appearance of the reconstruction changes: the spacing of the missing dimer trenches decreases, and the trenches become straighter. In this section, we investigate the effect on periodicity, while in the next section (Sec. 3.2) we model the straightening of the trenches.

The energy used for comparisons between different periodicities, $`E_p(\theta )`$, is $`10.35`$ eV/dimer for $`\theta =1`$ ML, $`10.08`$ eV/dimer for $`\theta =2`$ ML, and $`10.11`$ eV/dimer for $`\theta =3`$ ML. The small increase in this energy from 2ML to 3ML is surprising, but can be attributed to relaxation effects in the substrate (the third layer of Ge allows more relaxation than might be expected).

The results for different Ge ML coverage are given in table 1 and Figure 3. As shown in the figure, the optimal N value decreases when Ge coverage increases from 1ML to 2 ML, and again from 2 ML to 3 ML. This result indicates that increasing the Ge coverage results in more and more vacancies to release evolving stress and reduce strain energy. Figure 3 also shows that the increase in energy with N becomes steeper at the deeper coverages, indicating that the statistical distibution of N would become narrower with increasing Ge coverage (since the larger values of N become more expensive energetically), as also shown by STM. These results are broadly in agreement with what would be expected for an increasingly strained surface, and point towards the next stage in stress relief in this system: three dimensional islands (often known as “hut clusters”). The value of $`N=4`$ (or smaller) for 3ML is a little surprising — it is rather smaller than seen in experiment, or reported previously in modelling. The discrepancy with previous modelling is not surprising, as we use a different method of comparison with experiment. The discrepancy with experiment is somewhat more troubling, and will be discussed more in Section 4.

### 3.2 Ordering the $`(2\times N)`$ trenches

In experiments where the $`(2\times N)`$ reconstruction persists with increasing Ge coverage, as well as gradually reducing the periodicity N, the missing dimers become more ordered — i.e. they line up to form straight trenches. The obvious explanation for this is that it becomes energetically more favourable for the missing dimers to align with increasing Ge depth, and it is this phenomenon that we investigate here. Our strategy is to calculate the energy cost of disordering a trench of missing dimers as a function of the depth of Ge on the surface. We will introduce isolated kinks into straight trenches (i.e. a displacement of one or more dimers as we go from one dimer row to an adjacent one, as illustrated in Fig. 4) and calculate the energy difference between this configuration and a perfect trench.

We assume that only the top layer of Ge is defective, and use the optimal value of N for each system as calculated in section 3.1. We have performed kinking calculations on unit cells with two, four and six dimer rows, in order to find the interaction between adjacent kinks (i.e. in the two dimer row cell, the kinks will be as close as possible, while in the six dimer row cell, they will be separated by three dimer rows each). This effect has been found to be small.

The trench kinking energy of different depths of Ge coverage (1, 2 and 3ML) is shown in table 2. To give an idea of how likely kinks of different depths are, we have also calculated the distribution of the trench kink depths at 600C, assuming that these are distributed in thermal equilibrium. We have take the probability of a given kink depth, $`l`$, at any site to be given by a Boltzmann factor, $`P_le^{E_l/k_BT}`$, and have normalised the distribution to obtain a proportion of depths:

$$\underset{l}{}e^{E_l/k_BT}=1.$$
(2)
We can see that at all coverages, a straight trench (kink depth of zero) is most energetically favorable, becoming more so as the depth of Ge (and hence strain) increases. For a depth of 1ML of Ge, 45% of trench sites will kink, leading to the highly disordered appearance seen in STM. As the coverage increases, the probability of kinking decreases, so that at 3ML only about 6% of trench sites will kink by 2, and 72% of kink sites will remain straight. It is interesting to note that at all coverages the energy difference from a kink depth of 1 to 2 is larger than from a kink depth of 0 to 1, suggesting that there is a strong strain along the trench. Looking at Fig. 4, we can see this strain in the third layer silicon atoms (shown as dark circles) which are significantly more twisted in (c) than in (a). This is orthogonal to the strain that we have seen previously which was along the dimer row. These results emphasise the important role that strain plays in the morphologies of the $`(2\times N)`$ reconstructed surface in the process of growing Ge/Si(001).

### 3.3 The $`(M\times N)`$ reconstruction

Having considered the $`(2\times N)`$ reconstruction for multiple layers of Ge, we now move on to the two-dimensional $`(M\times N)`$ reconstruction (or the “patch” reconstruction), for 2 and 3 ML coverage of Ge, illustrated in Fig. 1 for 2 ML of Ge. For the case of 2 ML of Ge, we fixed N at 8 (from our previous calculation for a single monolayer) and changed the value M to find the most energetically favorable. We used the method described in Sec. 2.2 to calculate $`\zeta (M)`$ for the $`(8\times M)`$ reconstruction, with $`E_{\mathrm{p}_{\mathrm{strip}}}(\theta )=10.03`$ eV. These results are presented in table 3. We can see that the optimal M is 6, that is, of $`(8\times M)`$ reconstructions, the $`(8\times 6)`$ is most energetically favorable and most stable structure. When modelling the 3 Ge ML coverage case, we used a periodicity of 6 for N (again, considering our earlier results for the $`(2\times N)`$ surface with 2 Ge ML) and $`E_{\mathrm{p}_{\mathrm{strip}}}(\theta )=10.02`$ eV. As shown also in table 3, we found that the optimal M for $`(6\times M)`$ reconstruction is between 4 and 6. From these results we can see that the $`(M\times N)`$ reconstruction is still an effective strain relief mechanism, and that for three layers of Ge more relief is obtained (since the value of $`\zeta (M)`$ is more negative).

### 3.4 Comparing the $`(2\times N)`$ and $`(M\times N)`$ reconstructions

We have discussed in Sec. 2.2 how we will compare the stabilities of the $`(2\times N)`$ and $`(M\times N)`$ reconstructions relative to the perfect surface, and these results are given in table 4. We used values of $`E_\mathrm{p}(\theta )`$ and $`E_{\mathrm{p}_{\mathrm{strip}}}(\theta )`$ given above in Sections 3.1 and 3.3 to calculate these energies.

We see that the $`(2\times N)`$ reconstruction is more stable than the $`(M\times N)`$ for a 2ML coverage of Ge, and also for 3ML though by a much smaller margin. It is interesting that the stability of the $`(M\times N)`$ reconstructions for 3ML is increased when the second layer has a periodicity of $`N=8`$ (though the most stable periodicity for the $`(2\times N)`$ reconstruction with two layers of Ge is $`N=6`$). Nevertheless, these reconstructions are all stable relative to a perfect surface, and their formation under different growth conditions will depend on a subtle interplay of factors, such as: deposition rate; substrate temperature; and growth source.

## 4 Conclusions

Various important results can be drawn from our modelling of several layers of Ge on Si(001). First, we have shown that tight binding, and in particular $`𝒪(N)`$ tight binding, is sufficiently accurate to be used for modelling complex, strain-driven reconstructions such as those found in this system. Second, we have identified the trends behind the $`(2\times N)`$ reconstruction with increasing Ge depth, and shown that the trenches of missing dimers will tend to line up and straighten as the depth of Ge deposited increases. Third, we have explored the $`(M\times N)`$ “patch” reconstruction for different Ge coverages, contrasting it with the $`(2\times N)`$ reconstruction, and shown that, while it is less stable than the $`(2\times N)`$ reconstruction for all coverages of Ge considered, it becomes less unstable as the depth of Ge deposited increases.

Our results are based on two fundamental assumptions: full thermal equilibrium; and a lack of intermixing. We have already seen one area where our results appear not to be in exact agreement with experiment (in the value of N for the $`(2\times N)`$ reconstruction), and it is quite possible that these assumptions may be partly responsible for this discrepancy. The experimental observations generally take place on substrates which roughen as growth proceeds, leading to finite domains, which will have a considerable effect on the values of N observed. The values of M found for the $`(M\times N)`$ reconstruction show a similar trend towards values of $`M=4`$, though not as sharply (the value of $`M=6`$ for three layers of Ge is essentially degenerate with the value of $`M=4`$). At this size of reconstruction, some of the forces governing stability are due to the elastic constants of the system, which tight binding may not reproduce exactly. Nevertheless, our results reproduce and explain observed trends rather well, and may be considered as a limiting case, where the simplest assumptions apply: thermal equilibrium and complete Ge/Si segregation.

We thank the UK Engineering and Physical Sciences Research Council and the Royal Society (DRB) and the Chinese Government (KHL) for funding. We gratefully acknowledge useful conversations with Jaime Oviedo, Ilan Goldfarb and James Owen.


text:
# Origin of the photoemission final-state effects in Bi2Sr2CaCu2O8 by very-low-energy electron diffraction

## Abstract

Very-low-energy electron diffraction with a support of full-potential band calculations is used to achieve the energy positions, $`𝐊_{}`$ dispersions, lifetimes and Fourier compositions of the photoemission final states in Bi2Sr2CaCu2O8 at low excitation energies. Highly structured final states explain the dramatic matrix element effects in photoemission. Intense $`c`$(2x2) diffraction reveals a significant extrinsic contribution to the shadow Fermi surface. The final-state diffraction effects can be utilized to tune the photoemission experiment on specific valence states or Fermi surface replicas.

Photoemission (PE) spectra of Bi2Sr2CaCu2O8 (Bi2212) and other high-temperature superconductors can significantly deviate from the photohole spectral function $`A(𝐤,\omega )`$ due to an energy and $`𝐤`$ dependence of the matrix element (ME) . For example, the long debated dependence of the shape of the observed Fermi surface (FS) of Bi2212 on the photon energy $`h\nu `$ and experimental geometry has finally found its explanation in ME effects. A predictive use of these effects can help to resolve subtle details of the electronic structure of Bi2212 such as the bilayer splitting . Presently, analysis of ME effects employs computational simulations of the PE process . However, there remains an inherent uncertainty of a few eV due to the lack of any reliable knowledge about the energy location and lifetime of the PE final states. Additional complications are incommensurate (IC) structural modulations and antiferromagnetic (AFM) fluctuations, which manifest themselves in PE experiments, respectively, as diffraction replicas of the FS and so-called shadow FSs .

Very-Low-Energy Electron Diffraction (VLEED) gives direct access to the PE final states (see  and the references therein) based on the fact that, neglecting the electron-hole interaction within the sudden approximation, the PE final states are the time-reversed LEED states . In the VLEED elastic electron reflectivity spectra $`R(E)`$, the energies of the spectral structures give the final-state dispersions $`E(𝐤)`$, and their broadening the corresponding lifetimes $`\tau `$ expressed by the imaginary part $`Im\mathrm{\Sigma }`$ of the self-energy as $`\mathrm{}/\tau =2Im\mathrm{\Sigma }`$.

We here use VLEED to achieve first direct information on the PE final states in Bi2212 up to 45 eV above the Fermi level $`E_F`$ and elucidate its immediate implications for the PE experiment.

Experiment and band calculations – The VLEED experiment was carried out using a standard LEED setup operated in the retarding field mode . The total $`R(E)`$ integrated over all diffracted beams (plus structureless inelastic reflectivity) was measured in the target current $`T(E)`$. Intensities of individual beams were measured on the LEED screen using a CCD camera. Three samples with an oxygen content near optimal doping showed identical results.

The band calculations employed the full-potential APW+local orbitals method as implemented in the WIEN2k package . The generalized gradient approximation (GGA) was used to describe exchange-correlation within the density functional theory (DFT). The calculations employed the body-centered-tetragonal (bct) crystal structure and a basis of 2200 augmented plane waves. The obtained $`E(𝐤)`$, characterizing essentially the bulk states with infinite lifetime, was connected to the VLEED process based on Fourier expansion of the Bloch waves, which was used to evaluate the partial currents $`T_𝐤`$ into the crystal excited in each Bloch wave, and the integral $`T(E)`$ spectra . Finite electron lifetime was simulated by Lorentzian smoothing and damping of $`T(E)`$ .

The experimental structural parameters from  and  were tested in the calculations, with the latters yielding the best agreement with the VLEED experiment. While the valence bands changed insignificantly, the unoccupied $`E(𝐤)`$ was strongly sensitive to the structural parameters. Such a sensitivity suggests that any model of the final states, used in PE simulations, needs to be checked against the VLEED experiment.

Integral $`T(E)`$ spectra: Dispersions and lifetimes of the final states – The structures in $`T(E)`$ all reflect those in elastic $`R(E)`$, with the inelastic reflectivity giving only a featureless background. We start with analysis of the normal incidence spectrum in Fig.1(a). It corresponds, by parallel momentum conservation and assuming insignificant umklapp contributions from the IC superstructure, to the $`\mathrm{\Gamma }Z`$ direction of the Brillouin zone (BZ). A slope of $`T(E)`$ towards higher energies is due to an increase of the inelastic reflection. Its unusual steepness may indicate strong electron coupling in the cascade relaxation process. To emphasize the $`T(E)`$ structures, we subtracted a quadratic background, Fig.1(b).

The calculated $`T(E)`$ spectrum (minus a quadratic background) is also shown in Fig.1(b). A reasonable agreement with the experiment, especially considering the neglect of umklapp, confirms our theoretical model of the unoccupied states. The best averall agreement is achieved by shifting the calculations by $``$3 eV to higher energies. Besides the computational approximations, this can be attributed to excited-state self-energy corrections $`\mathrm{\Delta }\mathrm{\Sigma }`$ to the DFT-GGA exchange-correlation. The experiment suggests certain band dependence of $`\mathrm{\Delta }\mathrm{\Sigma }`$.

Guided by the calculated $`T(E)`$, we can trace the band structure origin of the experimental $`T(E)`$ structures. The $`T(E)`$ maxima are seen to reflect manifolds of so-called coupling bands (= dominant final states in PE) which are characterized by large partial $`T_𝐤`$ currents . The $`T(E)`$ minima reflect the gaps between such manifolds. Individual bands are however not resolved in $`T(E)`$ because their energy separation is much smaller than $`Im\mathrm{\Sigma }`$. Note that the coupling bands, apart from effective coupling to the incident plane wave, must have a 3-dimensional (3D) character with considerable $`k_{}`$ dispersion to transmit current into the crystal; they originate from the electronic orbitals oriented off-plane. The 2-dimensional (2D) bands, originating from the in-plane orbitals, are ineffective in both VLEED and PE.

Finite electron lifetime manifests itself in broadening and damping of the $`T(E)`$ structures. We estimated the energy dependence of $`Im\mathrm{\Sigma }`$ based on its approximate parametric relation to the dielectric function $`ϵ`$ as $`Im\mathrm{\Sigma }(E)=a+b_0^EIm\left[1/ϵ\left(𝐪=0,\omega \right)\right]d\omega `$ with $`ϵ\left(𝐪=0,\omega \right)`$ taken from the electron energy loss data . The parameters $`a`$ and $`b`$ were optimized to bring the calculated $`T(E)`$ to the best agreement with the experiment in broadening and damping of the spectral structures. The obtained $`Im\mathrm{\Sigma }(E)`$, which yields the calculated $`T(E)`$ in Fig.1 (b), is shown in the insert in Fig.1 (a) superimposed on $`ϵ\left(𝐪=0,\omega \right)`$. The increase of $`Im\mathrm{\Sigma }`$ results from two plasmon peaks in $`ϵ\left(𝐪=0,\omega \right)`$ associated with the CuO valence bands . This increase of $`Im\mathrm{\Sigma }`$ in the final states explains, for example, why the PE intensity ratio of the CuO bilayer bonding and antibonding states varies with $`h\nu `$ above 30 eV much smoother in the experiment than in simulations employing constant $`Im\mathrm{\Sigma }=`$1 eV  relevant closer to $`h\nu `$ of 20 eV.

Based on the experimental $`Im\mathrm{\Sigma }(E)`$ and calculated band dispersions, we found that the inelastic electron scattering limits the mean free path $`\lambda `$ around 3Å, decreasing with energy, in agreement with the electron energy loss data . In the final-state band gaps $`\lambda `$ is further reduced by elastic scattering off the crystal potential, resulting in a diffraction structure in $`\lambda (E)`$, particularly strong where $`Im\mathrm{\Sigma }`$ is small. Qualitatively, it follows the structure in $`T(E)`$ where the maxima reveal the 3D states enabling effective electron transport (accurate evaluation of $`\lambda `$ involves calculations of damped Bloch waves ). Predictive tuning $`\lambda `$ based on the experimental $`T(E)`$ can therefore be employed in the PE experiment to resolve signals from different atomic layers.

Experimental dispersion of the $`T(E)`$ spectra with the incidence $`𝐊_{}`$ varying in the $`\mathrm{\Gamma }X`$ azimuth (perpendicular to the IC superstructure) is shown in Fig.2 (left). The shading shows the energy intervals of the $`T(E)`$ maxima reflecting the coupling bands. This map gives therefore the surface-projected $`𝐊_{}`$ dispersions of the final states.

Our VLEED results give the first direct evidence of unoccupied quasiparticle dispersions in Bi2212, whereas the inverse PE studies near $`E_F`$ (see, e.g., ) found only dispersionless states. Furthermore, the highly structured non-parabolic dispersions illustrate a general fact that at low excitation energies the final states, especially for complex materials, dramatically deviate from a free-electron-like continuum due to multiple scattering on the crystal potential. This property of the final states is the essential origin of dramatic $`h\nu `$ dependences in PE.

The results for the $`\mathrm{\Gamma }Y`$ azimuth (along the IC superstructure) are shown in Fig.2 (right). By the fundamental lattice symmetry the $`\mathrm{\Gamma }Y`$ azimuth is equivalent to $`\mathrm{\Gamma }X`$, but this is complicated by the IC superstructure. The $`\mathrm{\Gamma }Y`$ dispersion image can then be viewed as the $`\mathrm{\Gamma }X`$ dispersion with its IC umklapp replicas superimposed and hybridized with the main dispersion. The observed image resembles therefore a smeared $`\mathrm{\Gamma }X`$ dispersion.

The dispersion in $`\mathrm{\Gamma }Y`$ comes at first glance as a surprise: Competing fundamental and IC potentials results in the collapse of strict periodicity (and the entire band structure concept) and thereby the parallel momentum conservation. As all states are now available for any incident $`𝐊_{}`$, one might expect that any dispersion in $`𝐊_{}`$ would be lost. However, the VLEED spectra depend, apart from bare availability of the states, on how their wavefunctions match the incident plane wave . Variation of the matching conditions with $`𝐊_{}`$ results in the observed dispersion. A similar appearance of competing periodicities in PE was elucidated in .

Our experimental $`T(E)`$ data yield thus the energy positions, $`𝐊_{}`$ dispersions and lifetimes of the PE final states. Apart from direct implications in the PE experiment such as tuning $`\lambda `$ in energy and $`𝐊_{}`$, this information can be used as a vital input in computational simulations of the ME effects instrumental for resolving details of the valence $`E(𝐤)`$ . ME changes are expected between the $`T(E)`$ maxima owing to a phase change in the final-state wavefunction through the band gaps.

Intensities of diffracted beams: hidden $`c`$(2x2) periodicity and final-state diffraction effects – A typical normal-incidence VLEED image is shown in Fig.3 (left). Due to the IC superstructure the diffracted intensity is arranged in lines along $`\mathrm{\Gamma }Y`$ , with non-vanishing intensity between the spots manifesting the fractal character of diffraction in presence of two competing periodicities. The intense spots, each with IC satellites, are seen to reflect, firstly, the reciprocal lattice of the fundamental lattice at $`(\pm 2\pi ,0)`$ and $`(0,\pm 2\pi )`$. Surprisingly, there are other intense spots $`(\pm \pi ,\pm \pi )`$ (arrows), although weaker in $`(\pi ,\pi )`$ and $`(\pi ,\pi )`$. Their symmetry reveals a hidden $`c`$(2x2) periodicity. This phenomenon is observed in electron diffraction for the first time.

$`I(V)`$ data are shown in Fig.3 (right) as the VLEED intensity profile along an IC diffraction line, as a function of energy and $`q_{(\pi ,\pi )}`$ momentum transfer. All diffraction spots show a significant and non-monotonous energy dependence beyond the kinematic picture. The $`c`$(2x2) spots are intense only below $``$20 eV, and in the normal LEED energy range they vanish.

To identify the origin of the $`c`$(2x2) periodicity, we first performed a kinematic simulation using the IC superstructure parameters from . It has however yielded only the fundamental-lattice spots with IC satellites. Beyond the kinematic picture, possible origins of the $`c`$(2x2) spots should be similar to the shadow FS with its $`c`$(2x2) periodicity: either two inequivalent Cu sites per plane in the full orthorhombic unit cell of Bi2212 , or short-range AFM correlations within the CuO planes  manifesting themselves in VLEED through the exchange scattering . However, the observed $`𝐊_{}`$ broadening of both $`c`$(2x2) spots and shadow FS seems too small compared to the short AFM correlation lengths found by neutron scattering. Finally, the $`c`$(2x2) spots can be peaks of multiple scattering on the competing fundamental and IC potentials. Whichever the origin, the diffraction process involves multiple scattering between the CuO planes, and the outer SrO and BiO planes with their IC modulations . The involvement of the CuO planes is corroborated the $`c`$(2x2) intensity reduction with energy which goes along with the decrease of $`\lambda `$.

Complementing the information on the final states obtained from $`T(E)`$, the VLEED intensity distribution reflects their Fourier composition: Being the time-reversed LEED state, the final state can be expanded over the surface reciprocal vectors $`\mathrm{\Phi }\left(𝐫\right)=_𝐠\varphi _𝐠\left(r_{}\right)e^{i\left(𝐤_{}+𝐠\right)𝐫_{}}`$, where the Fourier coefficients $`\varphi _𝐠\left(r_{}\right)`$ correspond, in vacuum, to the incident and diffracted plane waves and, in the crystal, to the Bloch wavefield. The VLEED intensities $`I_𝐠`$ are the vacuum asymptotics of $`\left|\varphi _𝐠\left(r_{}\right)\right|^2`$. As $`\varphi _𝐠\left(r_{}\right)`$ run across the surface continuously in amplitude and derivative, $`I_𝐠`$ reflect, approximately, the Fourier composition of the final state in the crystal.

In the PE experiment this information immediately characterizes the final-state diffraction effects. For example, in the FS mapping experiment a $`𝐠`$-diffraction replica of the main FS occurs whenever upon variation of $`𝐤_{}`$ the corresponing $`𝐤_{}+𝐠`$ Fourier component of the final state matches, in the parallel momentum, any Fourier component of the initial state, giving rise to non-zero overlap of the final- and initial-state wavefunctions. The PE signal is then proportional, apart from phases of the Fourier components and selection rules entering the ME , to the final-state $`\left|\varphi _𝐠\left(r_{}\right)\right|^2`$ and thus to its $`I_𝐠`$ asymptotics. Therefore, the PE signal of certain $`𝐠`$-diffraction replica follows, roughly, the $`𝐠`$-diffracted intensity observed in the VLEED image. The $`I(V)`$ data can then be used to tune the final-state energy to enhance or suppress certain diffraction replicas in the FS mapping experiment. Similar ideas can be exploited to single certain bands out of the valence band multitude based on their dominant Fourier component. Apart from Bi2212, such an approach can be applied to detail electronic structures of other complex materials, for example, perovskites.

Our VLEED data on the final-state diffraction effects leads to significant revisions in the current interpretation of some PE data on Bi2212. For example, the shadow FS with its $`c`$(2x2) periodicity is commonly considered as an intrinsic effect pertinent to the initial state. To check against the extrinsic effects due to the final-state diffraction, we have measured VLEED images with incidence $`𝐊_{}`$ set near the $`M`$ point where a large shadow FS signal is found in the PE experiment. The $`c`$(2x2) spots were found intense below $``$25 eV and gradually vanishing at higher energies. This indicates that at low $`h\nu `$, commonly used in the PE experiment, the shadow FS intensity has a large extrinsic contribution. This is corroborated by the absence of a hybridization gap between the main and shadow FS . The intrinsic shadow FS signal should be measured with higher $`h\nu `$ above 25 eV.

Finally, it should be noted that the final-state energies in PE, strictly speaking, differ from those observed in VLEED by the electron-hole interaction (excitonic) energy $`E_{eh}`$. This figure, characterizing deviations from the sudden approximation, can be determined by comparing characteristic points in the PE intensity behavior with the predictions of VLEED. For typical layered materials such an analysis has however yielded negligible $`E_{eh}`$ values, at least on the scale of 100 meV .

Conclusion – Using VLEED with a support of full-potential band calculations, we have achieved first direct information on the energy positions, $`𝐊_{}`$ dispersions, lifetimes and Fourier composition of the PE final states in Bi2212 at the excitation energies up to 45 eV above $`E_F`$. Implications of this information for the PE experiment extend from tuning $`\lambda `$ to predictive use of the ME effects to focus on specific valence states or FS replicas. In particular, intense final-state $`c`$(2x2) diffraction identifies significant extrinsic contribution to the shadow FS.

We thank M. Golden, A. Zakharov and E.E. Krasovskii for promoting discussions. The work was supported by Deutsche Forschungsgemeinschaft (grant Cl 124/5-1).


text:
# Blue Horizontal-Branch Stars and Simple Stellar Populations

## 1 Introduction

There have been many efforts to develop evolutionary population synthesis models (e.g. Bruzual & Charlot 1993; Worthey 1994; Buzzoni 1995; Vazdekis et al. 1996; Maraston 1998; Lee et al. 2000,2002) in order to analyze the integrated spectrophotometric quantities of globular clusters (GCs) and galaxies. Ultimately, the primary use of such software packages lies in attempting to derive the ages and metallicities of stellar populations such that quantitative insights into the processses governing galaxy formation can be uncovered. Age and metallicity drive the integrated spectrophotometry of stellar systems, and they also drive the resulting HB morphology of globular clusters. It is interesting therefore to pose the question as to whether this underlying morphology also manifests itself in the integrated spectrophotometry. There is reason to suspect that age is the global “second parameter” that controls HB morphology (after the “first parameter”, metallicity - e.g. Lee, Demarque & Zinn 1994; Sarajedini, Chaboyer & Demarque 1997; Rey et al. 2001; Salaris & Weiss 2002), although for some GCs a third (or more) parameter may be needed to explain their peculiar HB morphologies (such as the blue tail phenomenon - e.g. Recio-Blanco et al. 2002). In what follows, we describe the unique aspects of our recent models which allow for an accurate determination of the role played by blue HB stars in modifying111From that expected from the inclusion of the main sequence turnoff and red giant branch phases alone. the spectrophotometric indices of SSPs. The reader is directed to Lee et al. (2000,2002) for details.

## 2 Wavy Features in the Optical

We first construct a grid of simulated GCs spanning a range of HB morphologies \- in order to match the observed morphology of inner halo Galactic GCs (Galactocentric radius $``$ 8 kpc) at their currently favoured mean age ($``$12 Gyr) and age dispersion ($`\mathrm{\Delta }`$t=0 Gyr) a Reimers (1975) mass loss parameter $`\eta `$=0.65 was required. The Salaris et al. (1993) $`\alpha `$-element enhancement correction was applied to our models. The modified Gaussian mass distribution of Demarque et al. (2000) - with $`\sigma `$=0.02 M \- was used, shaping the widths of the “wavy” features seen in Figure 1. The temperature-sensitive B-V and H$`\beta `$ indices shown there illustrate the magnitude of the contribution of blue HB stars (thick wavy lines - c.f. Maraston & Thomas 2000) within $`\pm `$4 Gyr of system’s mean age. The resulting wavy features stem primarily from our exploration of parameter space 2$``$4 Gyrs older than the Milky Way GCs. The location of these wavy “peaks” is determined by the amount of mass loss adopted.

## 3 Usefulness of Far-UV Photometry

Perhaps the greatest value that the addition of far-UV photometry to pre-existing optical data provides is the ability to discriminate cleanly between young ($`<`$1 Gyr), intermediate-age (3$``$5 Gyr), and old ($`>`$14 Gyr), metal-rich ($``$0.8 $``$ \[Fe/H\] $``$ $``$0.4), clusters - this ability is demonstrated in Figure 2, where it can be seen that a significant far-UV to optical flux ratio is indicative of either very young or very old systems (for these metallicities). In contrast, intermediate-age clusters are relatively faint in the far-UV.

We suggest that the combination of far-UV and optical colours is a most promising chronometer for SSPs (with the admitted caveat that there are peculiar GCs such as NGC 6388 and 6441 with populations of extremely blue HB stars - e.g. Rich et al. 1997). If there are significant age differences between Milky Way and giant elliptical GC systems, however, in the sense of the latter being $`>`$2 Gyr older (as suggested by Lee et al. 2000), then we predict that the majority of metal-rich clusters will be far-UV bright compared to their Galactic counterparts. In this respect, the outcomes of the recent HST far-UV photometry of GC systems in Virgo ellipticals (PIDs#8643 & 8725) are highly anticipated; future work with Hubble’s Advanced Camera for Surveys will also be a fruitful avenue of research.222A program to which we would be more than happy to contribute … Finally, the upcoming GALEX mission will shed important light on the UV properties of the M31 globular cluster system, perhaps resolving the controversy regarding their supposed young, metal-rich, nature (e.g. Burstein et al. 1984; Barmby & Huchra 2000).

## 4 Nonlinear Colour-Metallicity Relations

Those Galactic GCs with minimal foreground extinction (filled circles: inner halo clusters; open circles: outer halo clusters) are contrasted with our models in the ($`B`$ $``$ $`V`$)o, ($`V`$ $``$ $`I`$)o, ($`M`$ $``$ $`T_1`$)o, and ($`C`$ $``$ $`T_1`$)o vs. \[Fe/H\] planes in Figure 3. The data are taken from Harris (1996, 1999 June version) and Harris & Canterna (1977). We believe this to be one of the first attempts at employing integrated broadband colours to derive relative age differences between subpopulations of Galactic GCs.

It appears that our models reproduce the differences between inner and outer halo clusters in the sense that the inner halo clusters are not only more tightly grouped along the isochrone than the more scattered outer halo counterparts, but also relatively older. This relative age difference between inner and outer halo Milky Way GCs is in agreement with that obtained by Salaris & Weiss (2002) from their homogeneous age dating of 55 clusters (as shown in Figure 4). The HB morphology is responsible for driving the non-linearity seen near \[Fe/H\] $``$ $``$1.6, demonstrating that the inclusion of morphological effects are important if one is to attempt subtle relative age-dating between SSPs.

As most of the linear colour-metallicity transformation relations (CMRs) in the literature are based upon the limited range in colour of Galactic GCs, caution must be employed when extrapolating their use to redder extragalactic systems (or if there are sizable age differences amongst globular clusters within and between galaxies). Recently, Kissler-Patig et al. (1998) provided a set of CMRs using the globulars of NGC 1399, but the apparent systematic offset has been questioned by Lee et al. (2002). It remains to be seen if the bimodal colour distributions of extragalactic GC systems can be interpreted as simple metallicity or age differences, or whether a more complicated interplay between metallicity and age must be invoked.

## Acknowledgements

H.-c. Lee sincerely thanks the organizers of this beautiful workshop for their hospitality. Support for this work was provided by the Postdoctoral Fellowship Program of the Korea Science & Engineering Foundation, the Creative Research Initiatives Program of the Korean Ministry of Science and Technology, and the Australian Research Council through its Linkage International (LX0346898) and Large Grant (A00105171) Programs.