text:
# Querying Databases of Annotated Speech

## 1 . Databases of Annotated Speech Recordings

Annotated corpora have been an essential component of research and development in language-related technologies for some years. Text corpora have been used for developing information retrieval and summarisation software (e.g. MUC , TREC ), automatic taggers and parsers and machine translation systems . In a similar way, annotated speech corpora have proliferated and have found uses across a rapidly expanding set of languages, disciplines and technologies \[www.ldc.upenn.edu/annotation/\]. Over the last 7 years, the Linguistic Data Consortium (LDC) has published over 150 text and speech databases \[www.ldc.upenn.edu/Catalog/\].

Typically, such databases are specified at the level of file formats. Linguistic content is annotated with a variety of tags, attributes and values, with a specified syntax and semantics. Tools are developed for each new format and linguistic domain on an ad hoc basis. These systems are akin to the databases of the 1960s. There is a physical representation along with a hand-crafted program offering a single view on the data. Recently, the authors have shown how the three-level architecture and the relational model can be applied to annotated speech databases . The goal of this paper is to illustrate our two approaches and to describe ongoing research on query algebras.

Before presenting the models we give an example of a collection of speech annotations. This illustrates the diversity of the physical formats and gives an idea of the challenge involved in providing a general-purpose logical characterisation of the data. The Boston University Radio Speech Corpus consists of 7 hours of radio news stories \[www.ldc.upenn.edu/Catalog/LDC96S36.html\]. The annotations include four types of information: orthographic transcripts, broad phonetic transcripts (including main word stress), and two kinds of prosodic annotation, all time-aligned to the digital audio files. The two kinds of prosodic annotation implement the system known as ToBI – Tones and Break Indices \[%\****␣arxiv.tex␣Line␣125␣\****www.ling.ohio-state.edu/phonetics/E\_ToBI/\]. We have added three further annotations: coreference annotation and named entity annotation in the style of MUC-7 \[www.muc.saic.com/proceedings/muc\_7\_toc.html\], and syntactic structures in the style of the Penn TreeBank . Fragments of the physical data are shown in Figure 1.

Coreference annotation (Figure 1, top left) associates a unique identifier to each noun phrase and a reference attribute which links each pronoun to its antecedent. The set of coreferring expressions is considered to be an equivalence class. Named-entity annotation (top centre) identifies and classifies numerical and name expressions. Penn Treebank annotation provides a syntactic parse of each sentence. The word-level annotation (bottom left) gives the end time of each word (a second offset into the associated signal data). The syllable annotation gives the Arpabet phonetic symbols (see \[www.ldc.upenn.edu/doc/timit/phoncode.doc\]). The tonal annotation provides time points and intonational units, and the part of speech annotation (bottom right) specifies the syntactic category of each word. This is but a small sample of the bazaar of data formats.

## 2 . Data Models for Speech Databases

Two database models for multi-layered speech annotations have been developed by the authors. The Emu model (Macquarie) organises the data primarily in terms of its hierarchical structure, while the annotation graph model (Penn) foregrounds the temporal structure. In separate work we demonstrate the expressive equivalence of the two models . Here we give a brief overview of both models. In the remainder of this paper we will consider mainly the annotation graph data model, while the Emu system serves as an example of a working speech database system.

### 2.1 . The Emu model

The Emu speech database system \[www.shlrc.mq.edu.au/emu\]  provides tools for creation, query and analysis of data from annotated speech databases. Emu is implemented as a core C++ library and a set of extensions to the Tcl scripting language which provide a set of basic operations on speech annotations. Emu provides a flexible annotation model into which a number of existing label file formats can be read.

The Emu annotation model is based on a set of *levels* which represent different types of linguistic data such as words, phonemes or pitch events. Each level contains a set of *tokens* which have one or more *labels* and optionally a start and end time relative to an associated speech signal. Within a level, tokens are stored as a partial order representing thier sequence in the annotation: each token may have zero or more previous and next tokens. The partial ordering must respect timing information if it is present in the tokens: that is, a token cannot follow a token with an later start time.

Within and between levels, tokens may be related by either *domination* or *association* relations. Domination relations relate a parent token to an ordered sequence of constituent child tokens and imply that the start and end times of the parent could be inferred from those of the children. Association relations have no in-built semantics and can be used for any application specific relation, such as that between a word and a tone target which denotes the point at which word stress is realised (Figure 2). Relations may be defined between any pair of levels which allows Emu to handle intersecting hierarchies such as that illustrated in Figure 2.

### 2.2 . The annotation graph model

A second general purpose model supporting multiple independent hierarchical transcriptions of the same signal data is known as the annotation graph . This model forms the heart of a joint initiative between LDC, NIST \[www.nist.gov\] and MITRE \[www.mitre.org\] to develop an architecture and tools for linguistic analysis systems (ATLAS), and an NSF-sponsored project between LDC, the Penn database group, and the CMU Psychology and Informedia departments, to develop a multimodal database of communicative interaction called Talkbank \[www.talkbank.org\].

Annotation graphs are labelled DAGs with time references on some of the nodes. Bird and Liberman have demonstrated that annotation graphs are sufficiently expressive to encompass the full range of current speech annotation practice. A simple example of an annotation graph is shown in Figure 3, for a corpus known as TIMIT . Annotation graphs (AGs) have the following structure. Let $`L=L_i`$ be the label data which occurs on the arcs of an AG. The nodes $`N`$ of an AG reference signal data by virtue of a function mapping nodes to time offsets $`T`$. AGs are now defined as follows:

###### Definition 1

An annotation graph $`G`$ over a label set $`L`$ and a timeline $`T`$ is a 3-tuple $`N,A,\tau `$ consisting of a node set $`N`$, a collection of arcs $`A`$ labelled with elements of $`L`$, and a time function $`\tau `$, which satisfies the following conditions:

1. $`N,A`$ is an acyclic digraph labeled with elements of $`L`$, and containing no nodes of degree zero;
2. $`\tau :NT`$, such that, for any path from node $`n_1`$ to $`n_2`$ in $`A`$, if $`\tau (n_1)`$ and $`\tau (n_2)`$ are defined, then $`\tau (n_1)\tau (n_2)`$;

Note that AGs may be disconnected or empty, and that they must not have orphan nodes. The AG corresponding to the Emu annotation structure in Figure 2, for the first five words of a TIMIT annotation, is given in Figure 3. The arc types are interpreted as follows: S – syntax; W – word; P – phoneme; T – tone; Imt – intermediate phrase; Itl – intonational phrase.

## 3 . Annotations as Relational Tables

Annotation data expressed in either the Emu or annotation graph data models can be trivially recast as a set of relational tables . For the purposes of this paper it is instructive to consider the relational form of annotation data in order to explore the requirements for a query language for these databases.

An annotation graph can be represented as a pair of tables, for the arc relation and time relations. The arc relation is a six-tuple containing an arc id, a source node id, a target node id, and three labels taken from the sets $`L_1,L_2,L_3`$ respectively. The choice of three label positions is somewhat arbitrary, but it seems to be both necessary and sufficient for the various annotation structures considered here.

We let $`L_1`$ be the set of types of transcript information (e.g. ‘word’, ‘syllable’, ‘phoneme’), and let $`L_2`$ be the substantive transcript element (e.g. particular words, phonetic symbols, and so on). We let $`L_3`$ be the names of equivalence classes, used here to model so-called ‘phonological association’. (This kind of association is discussed in depth in .) Let $`T`$ be the set of non-negative integers, the sample numbers. Figure 4 gives an example for the TIMIT data of Figures 2, 3.

We form the transitive closure of the (unlabeled) graph relation to define a structural (graph-wise) precedence relation using a datalog program:

```
s_prec(X,X) :-
s_prec(X,Y) :- arc(_,X,Y,_,_,_)
s_prec(X,Y) :- s_prec(X,Z),
               arc(_,Z,Y,_,_,_)

```

Now we further define a temporal precedence relation, where leq is the $``$ relation (minimally defined on the times used by the graph):

```
t_prec0(X,Y) :- time(X,T1),
                time(Y,T2),
                leq(X,Y)

t_prec(X,Y) :-  t_prec0(X,Y)
t_prec(X,Y) :-  t_prec(X,Z),
                t_prec0(Z,Y)

```

## 4 . Exploring Annotated Linguistic Databases

### 4.1 . General architecture

In our experience with the analysis of linguistic databases, we have found a recurrent pattern of use having three components which we will call query, report generation, and analysis.

The query system proper can be viewed as a function from annotation graphs to sets of subgraphs, i.e. those meeting some (perhaps complex) condition. The report generation phase is able to access these query results, but also the signals underlying the annotations. For example, the report generation phase can calculate such things as ‘mean F2 in signal S during time interval $`(t_1,t_2)`$.’ Each hit constitutes an ‘observation’ in the statistician’s sense, and we extract a vector of specified values for each observation, to be passed along to the analysis system. The analysis phase is then some general-purpose data crunching system such as Splus or Matlab.

This architecture saves us from having to incorporate all possible calculations over annotated signals into the query language. The report generation phase can perform such calculations, as well as compute properties of the annotation data itself. This seems to simplify the query system a good deal; now things like ‘count the number of syllables to the end of the current phrase’ (which we do need to be able to do) are tasks for the report generator, not the query system proper.

In general, the result of a query is a set of sub-graphs, each of which forms one matching instance. If we use the relational model proposed above, these would be returned as a result table having the same structure as the arc relation of Figure 4, but containing just the tuples which took part in each matching instance. We are then faced with the problem of how to differentiate the matching instances, for example, if we wished to collect together the word labels for the query ‘find all words dominated by noun phrases’ we need some way of treating each sub-graph separately. Hence, we would prefer the result to be a set of tables rather than a single table containing all matching tuples.

In a sense, then, the only role of the query is to define an iterator for the report generator over a set of sub-graphs of the overall annotation graph.

#### The Emu query language

The Emu query language uses simple conditions on token labels which match only tokens at a specified level, for example: Phonetic=A|I|O|U|E|V. These conditions can be combined by sequence, domination or association operators to constrain the relational structure of the tokens of interest. Examples of each are:

Find a sequence of vowel followed by stop at the phoneme level:
\[Phoneme=vowel -$`>`$ Phoneme=stop\]

Find Words not labelled x dominating vowel phonemes:
\[Word!=x ^ Phoneme=vowel\]

Find words associated with H\* tones:
\[Word!=x $`=>`$ Tone=H\*\]

The Word!=x query is intended to match any word in lieu of a query language construct which allows matching any label string.

Each query matches either a token or, in the case of the sequence query, a sequence of tokens. The result of a domination or association query is the result of the left hand side of the bracketed term; this can be changed by marking the right hand side term with a hash (#). Compound queries can be arbitrarily nested to specify complex constraints on tokens. As an example the following query finds sequences of stop and vowels dominated by strong syllables where the vowel is associated with an H\* tone target, the result is a list of the vowel labels with associated start and end times.

```
 [Syllable=S ^
     [Phoneme=stop ->
         [Phoneme=vowel => Tone=H*]]]

```

The result of an Emu query is a table with one entry per matching token:

```
database:timit
query:Phoneme!=x
type:segment
#
h#      0       147.5   fjsp0:sa1
sh      147.5   232.5   fjsp0:sa1
iy      232.5   325     fjsp0:sa1
hv      325     385     fjsp0:sa1
...

```
This table is used to extract any of the associated time-series data associated with the database, an operation usually carried out from an analysis environment such as Splus or XlispStat. Emu provides libraries of analysis functions for these environments which facilitate, for example, mapping signal processing operations over each token in a query result or overlaying plots of the time series data for each token.

Although this query system has proved useful and useable in the environment of acoustic phonetics research, it is now evident that there are a number of shortcomings which prevent it’s wider use. The query syntax is unable to express some queries, such as those involving disjunction or optional elements, and the query result is only really useful for data extraction. It is for these reasons that we are now looking more formally at the requirements for a query language for annotation data.

### 4.2 . A query language on annotation graphs

A high-level query language for annotation graphs, founded on an interval-based tense logic, is currently being developed and will be reported in a later version of this paper.

Here we describe a variety of useful queries on annotation graphs and formulate them as datalog programs. As we shall see, it turns out that datalog is insufficiently expressive for the range of queries we have in mind. Finding a more expressive yet tractable query language is the focus of ongoing research.

A number of simple operations, extending our two relations arc/6 and time/2, will be necessary for succinct queries. The first and most obvious is for hierarchy. Observe in Figure 3 that there is a notion of structural inclusion defined by the arcs. We formulate this as follows:

```
s_incl(I,J) :-
   arc(I,W,Z,_,_,_),
   arc(J,X,Y,_,_,_),
   s_prec(W,X), s_prec(Y,Z)

```

Now, since s\_prec is reflexive, so is s\_incl. Observe that nodes 3 and 6 in Figure 3 are connected by both an S/V arc and a W/had arc. The syntactic verb arc S/V should dominate the word arc W/had, but not vice versa. Therefore we need to have a hierarchy defined over the types. We achieve this with a (domain-specific) ordering on the type names:

```
type_hierarchy(word,syl)
type_hierarchy(syl,seg)

```

Now dominance is expressed by the predicate:

```
dom(I,J) :-
   arc(I,_,_,L1,_,_),
   arc(J,_,_,L2,_,_),
   type_hierarchy(L1,L2),
   s_incl(I,J)

```

In some cases it is necessary to have an intransitive dominance relation that is sensitive to phrase structure rules. For simplicity of presentation, we assume binary branching structures. The first of the rules below states that a sentence arc s will immediately and exhaustively dominate an np arc followed by a vp arc.

```
ps_rule(s,np,vp)
ps_rule(np,det,n)
ps_rule(vp,v,np)

```

Now we define immediate dominance over the syntax arcs syn as follows:

```
i_dom(I,J) :-
   arc(I,X,Z,syn,P,_),
   ps_rule(P,C1,C2),
   arc(J,X,Y,syn,C1,_),
   arc(_,Y,Z,syn,C2,_)

i_dom(I,J) :-
   arc(I,X,Z,syn,P,_),
   ps_rule(P,C1,C2),
   arc(_,X,Y,syn,C1,_),
   arc(J,Y,Z,syn,C2,_)

```

Another widely used relation between arcs is association. In the instance of the AG model in Figure 4, association amounts to sharing the value of $`L_3`$, as we saw in the tuples for dark and H\* in Figure 4. The assoc predicate simply does a join on the third label field:

```
assoc(I,J) :-
   arc(I,_,_,_,_,A),
   arc(J,_,_,_,_,A)

```

Finally, it is convenient to have a kleene star relation. Unfortunately in datalog we are unable to collect up the arbitrary length sequence it matches. Here we have it returning the two nodes which bound the sequence, which is often enough to uniquely identify the sequence in practice.

```
node(N) :- arc(_,N,_,_,_,_)
node(N) :- arc(_,_,N,_,_,_)

kleene1(X,X,_) :- node(X)
kleene1(X,Y,L) :-
   arc(_,X,Z,L,_,_),
   kleene1(Z,Y,L)

kleene2(X,X,_) :- node(X)
kleene2(X,Y,L) :-
   arc(_,X,Z,_,L,_),
   kleene2(Z,Y,L)

kleene3(X,X,_) :- node(X)
kleene3(X,Y,L) :-
   arc(_,X,Z,_,_,L),
   kleene3(Z,Y,L)

```

With this simple machinery we can start defining some annotation queries.

Find a sequence of vowel followed by stop at the phoneme level (assumes suitably defined vowel and stop unary relations):

```
vowel_stop(I,J) :-
   arc(I,_,Y,phoneme,V,_),
   arc(J,Y,_,phoneme,S,_),
   vowel(V), stop(S)

```

If we do not want both the vowel and the stop, but just the vowel, we could write:

```
vowel_stop(I) :-
   arc(I,_,Y,phoneme,V,_),
   arc(_,Y,_,phoneme,S,_),
   vowel(V), stop(S)

```

Find words dominating vowel phonemes:

```
strongWrdDomVowels(I) :-
    arc(I,_,_,word,s,_),
    arc(J,_,_,phoneme,V,_),
    vowel(V),
    dom(I,J)

```

Find words associated with H\* tones:

```
sylHtone(I) :-
    arc(I,_,_,word,_,A),
    arc(_,_,_,tone,h*,A)

```

Find stop-vowel sequences dominated by words in noun phrases where the word is associated with an H\* tone target.

```
stop_vowel_seq(I,J) :-
    arc(I,_,Y,phoneme,S,_), stop(S),
    arc(J,Y,_,phoneme,V,_), vowel(V),
    arc(W,_,_,word,_,_),
    arc(N,_,_,syn,np,_),
    dom(N,W), dom(W,I), dom(W,J),
    arc(T,_,tone,h*,_), assoc(W,T)

```

Find the intermediate phrase containing the main verb of a sentence:

```
imt_phrase(P) :-
    arc(K, _, _, syn, s, _),
    arc(J, _, _, syn, vp, _),
    arc(I, _, _, syn, v, _),
    i_dom(K,J), i_dom(J,I),
    dom(P, I),
    arc(P, _, _, imt, _, _)

```

Return the set of syllables between an H\* and an L% tone (inclusive).

```
syls(K) :-
    arc(_, _, N, tone, h*, A1),
    arc(_, N, _, tone, l%, A2),
    arc(I, _, N1, syl, _, A1),
    arc(J, N2, _, syl, _, A2),
    kleene1(N1, N2, syl),
    arc(K, N2, N3, syl,_,_),
    kleene1(N3, N4, syl)

```

The above query shows how the datalog model breaks down. We would like it to return sets of sets of syllable arcs. Instead it returns a flat set structure. In many cases we will know that some arc participating in the query expression can be used to recover the nested structure. For example, if the head of the above clause was changed from syls(K) to syls(I,K), then I will aggregate K in just the right way.

## 5 . Applying XML Query Languages to Annotations

It is worth briefly considering the suitability of existing XML query languages such as XML-QL  and XQL  for the domain of annotated speech. At first glance the problems we face querying annotated speech data are similar to those present with XML queries in that both present a hierarchical data model. A number of formulations of annotation data as XML are possible, indeed some projects make use of XML/SGML based formats entirely (e.g. MATE \[mate.nis.sdu.dklpq\], LACITO \[lacito.vjf.cnrs.fr/ARCHIVAG/ENGLISH.htm\]). XML can represent trees using properly nested tags, in the obvious way. In order to represent multiple independent hierarchies built on top of the same material one must construct trees using IDREF pointers. This idea was proposed by the Text Encoding Initiative  and recently adopted by the MATE project. We believe this approach is vastly more expressive than necessary for representing speech annotations, and we prefer a more constrained approach having desirable computational properties with respect to creation, validation and query.

The XQL proposal  describes a query language which is intended to select elements from within XML documents according to various criteria; for example, the query text/author returns all author elements that are children of text elements. The XQL data model ignores the order of elements within a parent element and has no obvious way to query for sequences of tokens.

The XML-QL proposal  provides for a data model where the order of elements is respected. A query for a word-internal vowel-stop sequence could be expressed as follows (assuming suitably tagged annotation data for TIMIT):

```
<word>
  <phoneme label=&vowel;/>
  <phoneme label=&stop;/>
</word>

```

The result of this query would have the following form:

```
<word label=had>
  <phoneme label=ae/>
  <phoneme label=dcl/>
</word>
<word label=dark>
  <phoneme label=ar/>
  <phoneme label=k/>
</word>

```

Queries which refer to two independent hierarchies, such as syntactic and intonational phrase structure, need to use joins. For example, to find words that are simultaneously at the end of both clauses and intermediate phrases, we could have the following query:

```
<intermediate>
  <word id=$i></>[end()]
</intermediate>

<clause>
  <word id=$i></>[end()]
</clause>

```

We assume the existence of some mechanism to pick out the last child element. The ID attribute ensures that the words are the same in each part of the join.

Perhaps either of these approaches could be made to work for a useful range of query needs. However they do not appear to be sufficiently general. For example, it is often useful to have query expressions involving kleene star: ‘select all pairs of consonants, ignoring any intervening vowels’ (CV\*C). Such queries may ignore hierarchical structure, finding sequences across (say) word boundaries. Using regular expressions over paths, XML-QL could provide access to strings of terminal symbols ignoring intervening levels of hierarchy. Yet it does not provide regular-expression matching over those sequences. Alternatively, sequences at each level of a hierarchy could be chained together using IDREF pointers, but it is unclear how we would manage closures over such pointer structures.

## 6 . Conclusions

Annotated speech corpora are an essential component of speech research, and the variety of formats in which they are distributed has become a barrier to their wider adoption. To address this issue, we have developed two data models for speech annotations which seem to be sufficiently expressive to encompass the full range of practice in this area. We have shown how the models can be stored in a simple relational format, and how many useful queries in this domain are first-order. However, existing query languages lack sufficient expressive power for the full range of queries we would like to be able to express, and we hope stimulate new research into the design of general purpose query languages for databases of annotated speech recordings.

## Acknowledgements

We are grateful to Peter Buneman, Mark Liberman and Gary Simons for helpful discussions concerning the research reported here.


format: text

text:
# Higher order variability properties of accreting black holes

## 1 Introduction

The continuum X-ray and $`\gamma `$-ray spectra of accreting black holes generally consist of a soft component that can be described as a multi-colored black body and a hard component that can be well fit by a cutoff power law with some additional flux due to Compton reflection (see e.g. Zdziarski et al. 1997; Gierlinski et al. 1997; Poutanen 1998). Thermal Comptonisation models (with an additional non-thermal component sometimes required) provide a physical basis for producing such a spectrum. A variety of models with differing geometries can fit the spectral data (see e.g. Poutanen 1998; Beloborodov 1999; Zdziarski 2000), so one must find a means of breaking the spectral degeneracies in order to differentiate between possible geometries of the flow and mechanisms for producing the high energy electrons that upscatter photons into the hard tail of the spectrum.

Another useful characteristic of accreting black holes, especially those in hard spectral states, is rapid variability (see e.g. Nolan et al. 1981; Negoro, Miyamoto & Kitamoto 1994; Vaughan & Nowak 1997). Several variability properties of their lightcurve have already been studied extensively, especially for Cygnus X-1, and place constraints on the underlying emission mechanism. For example, Fourier spectra and cross-spectra have been used to show that the spectral energy distributions of Cygnus X-1 and GX 339-4 (the two canonical persistent black hole candidates in the Milky Way) cannot be produced by Compton scattering in a static, uniformly dense cloud of hot electrons (Miyamoto et al. 1988), as had been previously suggested (Payne 1980). The improved photon statistics that are now possible with large area detectors like the Rossi X-ray Timing Explorer (RXTE) merit a re-examination of some of these prior analyses.

In a previous paper, we presented measurements of the cross-correlation function and autocorrelation functions of Cygnus X-1 at different energies which rule out all models in which the time lags between energy bands come from diffusion timescales through the corona. We also found that the autocorrelation and cross-correlation functions imply that the shape of individual shots in the hard state of Cygnus X-1 is such that the count rates rise more or less exponentially, then decay much more rapidly (Maccarone, Coppi & Poutanen 2000), the inverse of the FRED (fast-rise, exponential decay) behavior often cited in $`\gamma `$-ray burst studies (see e.g. Fenimore, Madras & Nayakshin 1996). This possibility can be evaluated and understood further by using the time skewness statistic, which measures the symmetry of a lightcurve as a function of time scale. Previous time skewness analyses showed that the light curve of Cygnus X-1 cannot be produced by purely exponential shots with either a rising or falling profile (Priedhorsky et al. 1979), but in our knowledge, this technique has not been applied since, despite immense advances in temporal resolution, photon statistics, and duration of observations. Here we present the results of time skewness analyses of RXTE data sets from Cygnus X-1 and GX 339-4 in their hard states. To complement this analysis, we also compute the Fourier bispectrum of the lightcurve data. The Fourier bispectrum (e.g. Mendel 1991) has been applied the spatial properties of the cosmic microwave background and large-scale structure (e.g. Cooray 2001) and to problems in other fields such as neurobiology (e.g. Bullock et al. 1995) and speech recognition (Fackrell 1996), but it has not previously been applied to astronomical time series data. With good quality data, however, we find that it provides another useful, higher order variability probe, determining whether the fluctuations on one time scale are coupled to the fluctuations on other timescales.

In section 2 we briefly describe the data set used and the data reduction and screening. In section 3, we define the statistics we use and present the observational results. In section 4, we discuss the mathematical form of the light curve implied by the observational results and consider the possible physical models that could produce the observed variability.

## 2 OBSERVATIONS

### 2.1 Reduction

Because of the need for large photon numbers to compute higher order variability statistics, we analyze the RXTE hard state data for the two brightest quasi-persistent low/hard state sources, Cyg X-1 and GX 339-4. The observations used are listed in Table 1. Light curves are extracted from the Proportional Counter Array (PCA) data in several different energy bands using the standard RXTE screening criteria of earth elevation greater than 10 degrees, pointing offset less than 0.01 degrees, having all 5 proportional counter units on, and having the standard amounts of time before and after South Atlantic Anomaly passages. In all cases, the time skewness is computed on a 0.125 second time scale. The photon statistics do not permit a computation on shorter timescales.

## 3 ANALYSIS

### 3.1 Time Skewness

#### 3.1.1 Definition

The time skewness measures the asymmetry of a lightcurve. We choose a slightly different definition than Priedhorsky et al. (1979) in order to non-dimensionalise the skewness and to make the normalisation independent of constant scaling factor changes to the count rate:

$$TS(\tau )=\frac{1}{\sigma ^3}\overline{[(s_i\overline{s})^2(s_{iu}\overline{s})(s_i\overline{s})(s_{iu}\overline{s})^2]},$$
(1)
where $`\tau `$ is defined to be $`u`$ times the time bin size, $`s_i`$ is count rate for the $`i`$th element of the light curve, and $`\sigma `$ is the standard deviation of the count rate. Essentially, the time skewness is a weighted average of the difference between two count rates measured at a given separation in time. If a light curve is composed of the summation of spikes, and the spikes typically rise more sharply than they fall on a particular time scale, the sign of the skewness will be negative; if the spikes fall more sharply than they rise, the skewness will be positive. We compute the time skewness for each observation separately, then present the weighted average of the individual measurements. The errors plotted are the rms variations of the individual measurements used in the averaging, and may be overestimated if there is intrinsic variation in the time skewness.

#### 3.1.2 Results

Throughout the rest of this paper the term “shot” will be used to describe a temporary increase in luminosity. In some cases in the literature this term has been used to describe models where the lightcurve consists of the random superposition of events with a single profile, duration, and intensity. Except when explicitly stated to have this meaning, we intend a more general definition of the word, and by referring to “shots”, we may imply a distribution of flare shapes and durations. Additionally, the term “rising shot” will be used to describe a shot whose duration is dominated by its rise time, and the term “falling shot” will be used to describe a shot whose duration is dominated by its decay time.

The 2-5 keV (channels 0-13 in the PCA) skewnesses from Cygnus X-1 and GX 339-4 in their hard states are presented in Figure 1. Both objects have a strongly positive skewness out to about 1 second and a strongly negative skewness between 1 second and the time scale on which their autocorrelation functions go to zero (5 seconds for Cygnus X-1 and 3 seconds for GX339-4). Beyond the autocorrelation function decay time scale, their skewnesses are consistent with zero. Whether these characteristics apply to all objects’ hard states is difficult to determine, because the other persistent hard state sources are much fainter than Cyg X-1 and GX 339-4. Furthermore, the transient hard state sources often show strong, relatively narrow QPOs in their hard states which complicate this analysis.

The skewnesses for other energy bands were also computed and are presented in Figure 2, with the error bars removed to make the figure more clear. In Cygnus X-1, the skewnesses are identical within the errors across energy bands. There seems to be some indication for a more symmetric overall lightcurve (i.e. a lower amplitude of the skewness) on long timescales in the higher energy bands, but it is not clear whether this is a real effect or due to the lower count rates at higher energy introducing more Poisson noise. In GX 339-4, the skewnesses become more symmetric as a function of increasing energy. This agrees well with the observation that the Fourier spectrum in GX 339-4 is consistent with being constant as a function of energy, while the Fourier time lags become longer for larger energy separations (e.g. Lin et al. 2000). We note that the large difference between the 15-45 keV lightcurve’s skewness and the lower energy bands’ skewnesses is due to the effects of Poisson noise artifically increasing the variance which is used to normalise the skewness values. This cannot be the case for the two lower energy bands, since the intermediate band has a count rate about 10% higher than the lowest energy band, and still has a lower skewness. Still the differences in the skewness as a function of energy are not particularly large, and it would be worthwhile to re-investigate this effect with the entire RXTE data set for GX 339-4 at the end of the mission.

The combination of the skewness measurements and the previous power spectrum results suggests that as the energy band increases, the rise timescale of a shot becomes slightly shorter while its decay timescale becomes somewhat longer. Thus the shots in GX 339-4 become more symmetric at higher energies and the normalisation of the skewness becomes smaller. Previous results found only that one-sided exponential shot models could be ruled out, and did not look to long enough timescales to see the sign change in the skewness (Priedhorsky et al. 1979). The fluctuations in the amplitude of the negative peak are of order the fluctutations in the skewness on 5-10 second timescales where the mean observed and expected values are $`0`$. We thus do not ascribe any physical meaning to the observed fluctations in the skewness within the peak.

The statistical properties of the time skewness are not well established. In order to ensure that our results are not due to meaaurement errors, we create simulated data sets of time-symmetric red noise with the phases at different frequencies uncorrelated with one another, using the algorithm of Timmer & König (1995), with a power spectral form similar to the observed power spectrum of Cygnus X-1 (a constant value below 0.1 Hz joined continuously to a $`f^1`$ power law from 0.1 to 1 Hz and then to a $`f^{1.7}`$ power law above 1 Hz, with an rms amplitude of 30%). One hundred lightcurves are simulated, each with a duration of 10 kiloseconds (shorter than the 30 kilosecond integration used in the Cygnus X-1 data set, which means that the real computation should be more robust than the simulated data computation). None of the simulated light curves show skewnesses on time scales of less than 1 second that could match the first peak (i.e. with absolute values greater than 0.02, the magnitude of the first peak). While three of the lightcurves show a skewness in the 2-4 second range with an absolute value of at least 0.04 (the size of the negative peak), there are no lightcurves with absolute values greater than 0.02 for the entire 2-4 second region, as is seen in Cygnus X-1. The mean value in the 2-4 second range is about 0.01, which are about the same size as the errors plotted in Figure 1. Thus, individual time bins from random light curves may occasionally (but infrequently) show skewnesses as large as the observed one, but they do not produce broad peaks in the skewness with depths as large as the observed depths. The randomly generated lightcurves do frequently have skewnesses half the size of those observed. This underscores the importance of long integrations for making skewness measurements, even when the photon statistics are good, since a short lightcurve, even with perfect photon statistics, can produce large skewness values despite having random Fourier phases at all frequencies.

The power spectrum for GX 339-4 is roughly the same as for Cygnus X-1, except that the high frequency cutoff occurs at about 3 Hz rather than 1 Hz. We simulate 40 kilosecond lightcurves (our observations are 46 kiloseconds here) with the power spectral shape for GX 339-4. As the observed skewness for GX 339-4 has a broad trough from 3 to 5 seconds at -0.015, we look for observations where the skewness has an absolute value of at least 0.015 over the entire 3 to 5 second range. We find no such lightcurves out of our 100 simulated lightcurves, and find only 3 lightcurves with any points with absolute values greater than 0.015 for any bin from 3 to 5 seconds. Thus again we find that while rarely, but occasionally, the depths of the troughs can be matched by the red noise lightcurves, troughs with the necessary depth and breadth cannot.

#### 3.1.3 Other spectral states

The time skewnesses were also computed for Cygnus X-1 in its transition state and in its soft state and for LMC X-1 which appears to have only a soft state. Since the rms variability amplitudes in the softer states are typically much lower than in the hard states (making variability harder to measure), and the characteristic variability time scales are shorter (making binning the data up to longer time scales problematic), the skewnesses computed were very noisy, and it was difficult to draw any useful conclusions from them. Thus we neither plot nor discuss these results.

### 3.2 Fourier bispectrum

#### 3.2.1 Definition

The bispectrum computed from a time series broken into $`K`$ segments is defined as:

$$B(k,l)=\frac{1}{K}\underset{i=0}{\overset{K1}{}}X_i(k)X_i(l)X_i^{}(k+l),$$
(2)
where $`X_i(k)`$ is the frequency $`k`$ component of the discrete Fourier transform of the $`i`$th time series (e.g. Mendel 1991; Fackrell 1996 and references within). It is a complex quantity that measures the magnitude and the phase of the correlation between the phases of a signal at different Fourier frequencies. Its value is unaffected by additive Gaussian noise, although its variance will increase for a noisy signal.

A related quantity, the bicoherence is the vector magnitude of the bispectrum, normalised to lie between 0 and 1. Defined analogously to the coherence function (e.g. Nowak & Vaughan 1996), it is the vector sum of a series of bispectrum measurements divided by the sum of the magnitudes of the individual measurements. If the biphase (the phase of the bispectrum) remains constant over time, then the bicoherence will have a value of unity, while if the phase is random, then the bicoherence will approach zero in the limit of an infinite number of measurements. The bicoherence is defined as:

$$b^2(k,l)=\frac{\left|X_i(k)X_i(l)X_i^{}(k+l)\right|^2}{\left|X_i(k)X_i(l)\right|^2\left|X_i(k+l)\right|^2}.$$
(3)
This analysis will concentrate on the bicoherence, since statistically significant measurements of it are easier to make. The denominator of this expression has a dependence on the amount of Gaussian noise, so a correction will have to be made when the signal is noisy. Additionally, the bicoherence, being bound between 0 and 1, has a non-zero mean value due to errors even when the phases of the Fourier transforms are uncorrelated. The mean value and the standard deviation of a bicoherence estimate based on $`K`$ distinct Fourier transforms are both $`\frac{1}{K}`$. The properties of the bispectrum and the bicoherence as they relate to time series analysis are well-reviewed by Fackrell (1996).

#### 3.2.2 Computation method and results

The lightcurves are split into segments of 4096 elements (since the time resolution of the lightcurves varies, the lengths of these segments in absolute amounts of time will also vary). Fourier transforms are computed for each element. The results are then combined according to equation (2) to estimate the bicoherence. The data used to compute the bispectra are listed in Table 2.

Since the bicoherence is a function of two independent variables, one could presumably gain the most information from a 3 dimensional plot, or a two dimensional contour plot. Since we find from inspection that, for these data, no individual peaks stand out in such plots, the data are rebinned to present the results only as a function of the sum of the two lower Fourier frequencies. This allows for a simple two-dimensional plot, allows for significant data rebinning to increase the statistical significance of the results, and does not remove any obvious underlying trends. Similar techniques are frequently used in analyses of spatial bispectra of the microwave background, although in that case the computations are done in terms of spherical harmonics rather than Fourier frequencies. We caution the reader who may wish to try this technique on other time series that for systems with strong, relatively narrow quasi-periodic oscillations, rebinning the frequencies to one dimension may result in the loss of important signals. For example, we have examined a few observations of GRS 1915+105 where the 1-10 Hz QPOs show harmonic structure and have found that the bicoherence can be quite strong for the pairs of frequencies that add up to the frequencies of the higher harmonics, but that this effect can be seen clearly only from a three dimensional plot.

We have plotted the results for the bicoherence measurements for Cygnus X-1 (from 2-5 keV; channels 0 to 13) and GX 339-4 (from 2 to 8 keV; channels 0 to 25) in Figures 3 (a) and (b) respectively. In all cases, the bispectrum fits a broken power law with a negative index for high frequencies; for lower frequencies, it has a nearly constant value. The non-zero bicoherence indicates that the variability at low frequencies is coupled to the noise at high frequencies, providing an additional piece of strong evidence that a single shot model cannot match the data. Furthermore the phases of the Fourier components are more strongly correlated at low frequency than they are at high frequency.

We attempted similar analyses for data from accreting neutron stars, but the photon statistics were not sufficient to find interesting results. There is not sufficient RXTE data on any single neutron star in the island state of an atoll source (the analog of the hard state for black hole sources) to make use of this technique. Likewise, LMC X-1 (the most variable and best observed soft state source) lacks sufficient signal to noise to make use of the bispectrum. For objects with strong quasi-periodic oscillations, the analysis is significantly more complicated, as some frequencies may represent the harmonics of periodicities found at lower frequencies. Such analysis is outside the scope of this paper which aims merely to characterise random noise components of the power spectra of accreting objects.

## 4 Discussion

We now consider possible models and attempt to determine how they match the key observational results - that (1) the skewness is positive on short timescales, (2) the skewness changes sign and becomes negative on timescales of about 1 to 5 seconds before going to zero on long timescales and (3) that the bicoherence has a “shelf” at low bifrequency. We also will make use of the result of Priedhorsky et al. (1979), reproduced here, that the normalisation of the skewness found in the data is much smaller than the normalisation of the skewness found from a one-sided exponential shot distribution. An exact reproduction of the time skewness is likely to be a quite complicated problem numerically, as the shot shapes and arrival times may take non-standard functional forms. We first discuss toy models that approximate the mathematical form of the skewness then discuss recent theoretical models which might be expected to produce this mathematical form. Finally, we show that the bicoherence results agree with the general mathematical form used to match the skewness, providing an independent confirmation of our picture. In the models we consider, we invoke the “shot noise” hypothesis - that the variability consists of discrete events. Our conclusions that the lightcurve shows a different sense of asymmetry on different timescales are independent of this assumption. The specific interpretations below are not so clearly robust.

One possible means of producing a skewness function which has both positive and negative signs is to create a light curve which has both rising and falling shots, with different time scales for the rising shots than for the falling ones. This would seem to imply two physical mechanisms for producing shots. Furthermore, to produce hard lags on all time scales from this distribution, as is required from the Fourier cross spectrum, one would need the rising shots to be shorter with increasing energy while the falling shots would have to be longer with increasing energy. This would produce an autocorrelation function that becomes wider with increasing energy on long time scales, in violation of observational constraints. We were unable to find a distribution of shot intensities and time scales for this scenario that reproduced the observed skewness.

A second model one might consider would be to have a random superposition of identical two-sided shots, with a rise time different from the fall time. Two sided shots, however, fail to reproduce to observed skewness measurements. They are always dominated by the longer timescale, and actually fail to produce skewnesses which have both substantial positive and negative portions on different time scales.

A third possibility is to have all the shots be rising exponentials with sharp cutoffs (inverse FREDs), while the probability of having a shot occur at a given time is given by an exponentially falling function (an actual FRED, i.e., the mean shot rate is proportional to $`\lambda \mathrm{exp}[(tt_0)/\tau ]`$ for $`(t>t_0)`$. Thus the rising shots are embedded in a falling envelope. In this paper, the envelopes will be modeled by drawing the arrival times of the peaks of the shots from an exponential deviate; more complicated envelope functions are possible but are not considered here. Figure 4(a) illustrates the difference between rising and falling shots, while Figure 4(b) shows rising shots embedded in a falling envelope.

An “envelope” model with one-sided exponential rising shots with a rise time of 0.5 seconds and a one-sided exponential envelope with a decay time of 4 seconds matches the qualitative trends of the observed skewness, but the overall normalisations of the skewness are a factor of $`10`$ too large, indicating that this model assumes much more asymmetry than is present in nature. We thus modify the model by giving both the envelope and the shots two-sided profiles. The resulting model lightcurve has a skewness that matches the observations within factors of a few. For the time scales on which the individual shots occur, the skewness is positive. For time lags of duration between the individual shot time scale and the “envelope” time scale, the skewness is negative. For time lags longer than the envelope scale, the skewness is zero. Skewness computations for a simulation with shots having an 0.5 second rise time scale and an 0.1 second decay time scale embedded in an envelope with a 2.7 second rise time scale, a 3.0 second decay time scale and 53% of the shots in the decaying portion (so that the rising and decaying portions of the envelope will form a continuous function) are computed. Each model lightcurve contains 25 shots. Since the process is random, 149 lightcurves have been computed and the results have been averaged to improve the statistics. The simulated skewness for this model is plotted in Figure 5 (a) and a sample simulated lightcurve is plotted in Figure 5 (b). We have also computed the skewness of a simulated lightcurve made of a random superposition of the 25 shot lightcurves. We find that the smearing of the shots induced by this superposition does not affect the results substantially. In order to fit the skewnesses more exactly, one needs to modify the parameters of the simulation by varying the ratios of rise time to decay time. A model in which the shot and envelope time scales are allowed to vary with some small range around an average produces essentially the same skewness results as the single time scale models, and if the number of shots per unit time is too high, then the individual shots blend together so much that the skewness reflects only the envelope. The idea that the shots cannot be independent stochastic events was also suggested by Uttley & McHardy (2001) to explain the constancy of RMS variability as a function of count rate in Cyg X-1, SAX J1808.4-3658, and several Seyfert galaxies.

A final possibility is related to the third one. One can have the shot arrival times be independent of one another, while their intensities are correlated with one another. If large shots are followed only by smaller shots, then the general form of the skewness also can be reproduced.

Such a mathematical form can be reproduced by two recent models - the self-organised criticality (SOC) model (Takeuchi, Mineshige & Negoro, 1995; Takeuchi & Mineshige 1997) and the pulse avalanche model (Poutanen & Fabian 1999). The skewness of a simulated light curve that borrows the major results of the SOC model is presented in Figure 6a, and a segment of the simulated light curve is plotted in Figure 6b. In this simulation, an energy reservoir is posited and shot arrival times are presumed to be random. Each shot results in the release of a random fraction of the energy reservoir. The duration of the shot is related to the energy released at $`\tau E^{0.4}`$. The energy reservoir is refilled at a constant rate. Thus when a large shot occurs, the reservoir’s energy level will be depleted and subsequent shots will be smaller until the reservoir is re-filled. The shots are given Lorentzian profiles and their $`\tau `$ values are given by the ratio of the reservoir’s energy level to the mean reservoir energy level times a random value drawn from a uniform distribution between 0.05 and 0.15 seconds. A 1000 second simulated light curve is produced, allowing us to average over many reservoir-filling timescales.

Current applications of the SOC model to black hole binary light curves fail to match several well established observational constraints for Cygnus X-1 and GX 339-4. They assume a shot with the reverse time sense of those seen in the observations. This is not an essential component of the model since the shot profile is put in by hand to match the power spectrum, and our simulation has already changed this property of the shots and assumed the shots decay rapidly rather than rise rapidly. Additionally, this model does not attempt to explain the time lags seen in black hole binaries in the hard state. Both these problems could probably be solved, at least qualitatively, by assuming that the gas in a shot event gets hotter as it falls in, then abruptly stop emitting light when it reaches the event horizon of the black hole. This would result in a shot with a rising profile, a sharp turnoff, and hard lags.

More troubling is the inability of the ADAF/SOC model to match the observed luminosity. Because it assumes the standard parameters from Narayan & Yi (1995), the total amount of energy in electrons in the corona is small (of order $`10^{34}`$ to $`10^{35}`$ ergs), the accretion rate must be small, and the accretion flow must radiate inefficiently. This makes it difficult to supply the required luminosity for Cygnus X-1 and GX 339-4 in their hard states, since the hard states of these sources are more luminous than those of most other accreting black hole candidates. Thus the models of Takeuchi & Mineshige (1998) predict a luminosity of only $`10^{34}`$ ergs/sec, while the luminosity of Cygnus X-1 in its hard state is about a factor of 1000 higher than that value.

Still, the skewness of the lightcurve is fairly well matched by the SOC model. This success may suggest that some conceptually similar process - an energy reservoir of some sort drives the variability; when that reservoir is depleted, production of large shots is suppressed. We also note that a gradual mass flow in addition to the self-organized critical flow is required to produce an $`f^1`$ power spectrum rather than a steeper $`f^2`$ decay. Whether some form of the SOC model could reproduce the lightcurve of Cygnus X-1 given a set of coronal parameters which reproduces the spectrum more accurately (such as that of Esin et al. 1998) remains to be seen. At the present, the problem is with the specific parameters used in the ADAF realisation of the SOC model and not with the SOC model itself. We note that our simulated “envelope” model has a similar problem of having its power spectrum fall off too rapidly at high frequencies. To fix the problem, a sharper than exponential envelope function would be needed.

The magnetic flaring model of Poutanen & Fabian (1999) is an alternative mechanism for producing the both the hard emission and the variability seen in low/hard state objects. This model produces the time lags through temperature changes in an active region. Each shot has a probability of causing another shot to occur and the mean delay time before the emission of the stimulated shot is proportional to the energy of the stimulating shot. The requirement that the magnetic energy reservoir be refilled before another shot occurs could be the cause of this delay. The possibility that a magnetic energy reservoir is required to heat the corona in accreting systems has been discussed recently (Merloni & Fabian 2001).

Another recent observational result that can perhaps be explained in terms of these multiple time scales, either due to reservoir effects or exponential envelopes of shots, is the decomposition of the Fourier power spectrum into multiple Lorentzian peaks (e.g. Nowak 2000). This approach typically yields between 3 and 7 peaks in the power spectrum, but the basic idea that the variability has multiple characteristic time scales is present in both our results and the “many-QPO” interpretations of the power spectra. The strongest peaks in the Fourier spectrum may correspond to the time scales presented here and the other peaks in the Fourier spectrum may be unresolvable with current skewness measurements or may require yet higher order moments of the light curve to be computed.

Finally, assuming that a shot-noise type model is correct, it is noted that in Cygnus X-1, these results present evidence for a significant constant (or quasi-constant) flux component. Since the shots have now been fairly clearly demonstrated to be shorter at higher energies, the rms variability expected due to the shots should rise as a function of energy. But in the observations (e.g. Nowak et al. 1999), the rms decreases as a function of increasing energy. In order to reconcile these two observations, one must include an additional component with a significantly harder spectrum than that emitted by the shots, or we must assume that the shots have different properties at different energies (e.g. different relative amplitudes or different arrival frequencies). In the context of a shot model, the variable relative amplitude hypothesis would specifically require that the brightest shots would have softer spectra than the faintest ones, while the variable arrival frequency would require that some shots show up in the hard energy bands without contributing to the low energy bands. Given the observation that the coherence function is greater than 0.95 on all timescales from 0.01 Hz to 10 Hz, but drops significantly above 10 Hz (Nowak et al. 1999), the hypothesis that there are more shots at high energies than at low energies can be restricted to the case where the extra shots at high energies contribute almost no power on timescales longer than 0.1 seconds. While the near unity cross-correlation functions can be used to constrain how much the shots’ spectra can vary (Maccarone, Coppi & Poutanen 2000), rather large differences in shot spectra can be allowed as long as the peaks are aligned. Positing a constant component with a harder spectrum than the variable component would seem to contradict the results of Uttley & McHardy (2001) who found evidence for a constant component with nearly the same spectrum as the variable component. Thus it seems that the picture where the brightest shots have the softest spectrum is the most likely one. The possibility that the reflection component contributes to the drop in rms amplitude at the highest energies (those $`>`$ 13 keV, which were not studied by Uttley & McHardy) remains viable.

Essentially the same implications are found from the bicoherence measurements as from the time skewness measurements. We find that the $`f^1`$ dependence of the bicoherence function at high frequencies cannot be well matched by a single shot model. In Figure 7 we show the results from a model bicoherence computation for a simulated lightcurve with an envelope that rises instantaneously and decays exponentially on a timescale of four seconds, which is composed of shots that rise on an exponential timescale and decay instantaneously. The rise times of the shots are distributed uniformly between 0.5 and 1.5 seconds. The resulting bicoherence measurements are plotted in Figure 7, and match the observations reasonably well. Because of the computational intensity of bicoherence calculations and the uncertainties about whether the shot shapes should be exponentials, stretched exponentials, or some other more complicated profile, we cannot find an exact set of parameters for the shots and their envelopes.

## 5 Summary

We have shown that the hard states of Cygnus X-1 and GX 339-4 have qualitatively similar skewness and bispectrum properties. They cannot be matched by any of the simplest single shot models. We suggest the possibility that either a model with a falling envelope of rising flares or a model with an energy reservoir that suppresses large shots following other large shots can reproduce the data, and that only models in which the arrival times and/or luminosities of the shots are correlated with one another can reproduce the skewness function which has both positive and negative parts. Whether the reservoir is one of accreting mass or of magnetic energy cannot be determined with the current data. Of the current variability models, the self organised criticality models or the magnetic flaring pulse avalanche models seem most likely to be able to fit this data without major modifications.

## 6 Acknowledgments

We thank Rick Rothschild for pointing out that the time skewness could be used to test whether shots are symmetric and Roger Blandford for suggesting the use of the bispectrum to augment the skewness results. We are grateful to Phil Uttley for useful discussions regarding the constant component in Cygnus X-1. We thank the anonymous referee for making several useful suggestions. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center.


text:
# Heterogeneous aggregation in binary colloidal alloys

## I Introduction

Soft condensed matter are materials which can easily be deformed by external stresses, electric and magnetic fields, or even by thermal fluctuations . These materials typically possess structures which are much larger than atomic or molecular scales; the structure and dynamics at the mesoscopic scales determine the macroscopic physical properties. The goal of our research is to model and understand this relationship. These materials can be synthesized by means of colloidal self-assembly. In what follows, we will describe a model colloidal system in which the particles interact through type-dependent forces, so as to study the relationship between the interaction at the mesoscopic scales and the macroscopic properties.

There is a long-standing yet fundamental question in physics: Given all the interactions between the particles, what will be the resulting structure that the particles will form? For monodisperse case, usually close-packed structures will be formed, e.g., hexagonal structure in two dimensions and face-centered cubic or hexagonal close-packed structures in three dimensions. By allowing two different types of particles, driven by two types of forces, we are prepared to show that a rich variety of possible structures will be formed, as we will study in more detail below.

The plan of the paper is as follows. In the next section, we will discuss the molecular dynamics simulation method. In section III, we will present the results. Discussion and conclusion on our results will be given.

## II Molecular Dynamics Simulation

The colloidal system studied in this work is a two-dimensional one consisting of circular particles of the same diameter $`d`$, suspended in a viscous fluid. There is a pairwise interparticle force $`𝐅_{ij}=U(𝐫_{ij})`$ between particles $`i`$ and $`j`$, where $`U`$ is a potential and $`𝐫_{ij}=𝐫_i𝐫_j`$. The equation of motion of a particle is given by

$$m\frac{d^2𝐫_i}{dt^2}=\underset{ji}{}𝐅_{ij}w\frac{d𝐫_i}{dt},$$
(1)
where the first term on the right-hand side describes the total force acting on the particle $`i`$, while the second term denotes the viscous drag (with coefficient $`w`$), exerted by the fluid. We have neglected the Brownian motion, which is a valid assumption for mesoscale objects. The variables can be rescaled as $`t=t_0t^{}`$, $`𝐫=d𝐫^{}`$, and $`𝐅=F_0𝐅^{}`$, with $`t_0=m/w`$ and $`F_0`$ being the typical magnitude of the interparticle force. The rescaled equation of motion can be written as

$$\ddot{𝐫}_i^{}=A\underset{ji}{}𝐅_{ij}^{}\dot{𝐫}_i^{},$$
(2)
where $`A=F_0t_0/w`$. In a highly viscous medium, $`A`$ is very small and the particles are in overdamped motion, i.e., $`\ddot{𝐫}_i^{}0`$. In this case, Eq. (2) reduces to

$$\dot{𝐫}_i^{}=A\underset{ji}{}𝐅_{ij}^{},$$
(3)
We further define $`\stackrel{~}{t}^{}=At^{}`$ and obtain the reduced equation of motion for our simulation:

$$\frac{d𝐫_i^{}}{d\stackrel{~}{t}^{}}=\underset{ji}{}𝐅_{ij}^{}.$$
(4)
Thus we have a set of simultaneous first-order differential equations to solve. In the overdamped situation, the final structure is independent of $`A`$, but the time for the structure formation is inversely proportional to $`A`$. The initial configurations consist of $`N`$ particles randomly dispersed in a square simulation cell, with periodic boundary conditions being imposed in two directions. Thus, we include the forces acting on a particle due to all other particles in the simulation cell as well as their periodic images, calculated within a square centered at that particle which is of the same size of the simulation cell. The simulations are run for a sufficiently long time so that the particles have no further movements as time goes on.

The simulation cell is set to the scales $`25d\times 25d`$ and we ran simulations of area fraction of particles $`\mathrm{\Phi }=0.13`$, $`0.19`$, $`0.25`$, $`0.31`$, $`0.37`$ and $`0.44`$ corresponding to the number of particles $`N=100`$, 150, 200, 250, 300, 350 respectively. The simple Euler algorithm has been used throughout the simulation with a time step $`\delta \stackrel{~}{t}^{}=5\times 10^4`$. Larger time steps may result in unphysically large velocities because of serious overlapping between the particles.

We consider two cases. The first case consists of one type of particles, interact through attractive force, while the second case consists of equal number of different types of particles A and B. The initial positions and velocities of the particles are imposed by a random distribution inside the simulation cell without particle overlapping. It is difficult to generate for $`N`$ larger than 350 because of serious particle overlapping.

### A One type of particles with attractive force

In this case all particles attract each other at short distance and repel each other at long distance. The potential between the partilces is

$$U(r)=u_0\mathrm{exp}((rd)/\xi _0)u_1\mathrm{exp}((rd)/\xi _1).$$
(5)
The first term denotes a repulsion, with a range $`\xi _0`$ and magnitude $`u_0`$. The second term denotes an attraction, with a range $`\xi _1`$ and magnitude $`u_1`$. We set $`\xi _0>\xi _1`$ and $`u_1>u_0`$ so that the particles attract each other at short distance. In order to avoid too much overlap when the particles approach and finally touch, we add a steep repulsive potential of the form $`b/r^{11}`$ to the interparticle potential. Thus, the force is calculated by the negative gradient of the potential

$$F(r)=\frac{11b}{r^{12}}+\frac{u^{}}{\xi _0}\mathrm{exp}((rd)/\xi _0)\frac{u_1}{\xi _1}\mathrm{exp}((rd)/\xi _1)$$
(6)
where

$$u^{}=u_0\frac{b}{d^{11}},\mathrm{and}b=\frac{u_0/\xi _0u_1/\xi _1}{1/\xi _0d^{11}11/d^{12}}.$$

In this way, the two particles attract each other when they approach, while the force vanish when $`r=d`$, and repel each other when they overlap. If we rescale the force by $`b=F_0d^{12}b^{}`$, $`u^{}=F_0du^{}`$, $`u_0=F_0du_0^{}`$ and $`u_1=F_0du_1^{}`$, we obtain the dimensionless force:

$$F^{}(r)=\frac{11b^{}}{r_{}^{}{}_{}{}^{12}}+\frac{u^{}}{\xi _0^{}}\mathrm{exp}((r^{}1)/\xi _0^{})\frac{u_1^{}}{\xi _1^{}}\mathrm{exp}((r^{}1)/\xi _1^{}),$$
(7)
where

$$b^{}=\frac{u_0^{}/\xi _0^{}u_1^{}/\xi _1^{}}{1/\xi _0^{}11}\mathrm{and}u^{}=u_0^{}b^{}.$$

Here $`\xi _0^{}`$, $`\xi _1^{}`$, $`u_0^{}`$ and $`u_1^{}`$ are adjustable parameters. In this simulation, the parameters used are $`\xi _0^{}=0.5`$, $`\xi _1^{}=0.3`$, $`u_0^{}=5`$, and $`u_1^{}=6`$.

### B Two type of particles

The system consist equal number of A and B particles of total area fraction $`\mathrm{\Phi }`$. For unlike particles, there are short range attractive and repulsive forces between the particles, same as for one type of particles above. The potential between two particles is

$$U_{}(r)=u_0\mathrm{exp}((rd)/\xi _0)u_1\mathrm{exp}((rd)/\xi _1),$$
(8)
while for the same type of particles, there is always a repulsive potential between them

$$U_+(r)=u_0\mathrm{exp}((rd)/\xi _0).$$
(9)
Again, we impose a steep repulsive potential of the form $`b/r^{11}`$ to the interparticle potential. Likewise, the rescaled dimensionless force are:

$$F_{}^{}(r)=\frac{11b^{}}{r_{}^{}{}_{}{}^{12}}+\frac{u^{}}{\xi _0^{}}\mathrm{exp}((r^{}1)/\xi _0^{})\frac{u_1^{}}{\xi _1^{}}\mathrm{exp}((r^{}1)/\xi _1^{})$$
(10)
$$F_+^{}(r)=\frac{11b^{}}{r_{}^{}{}_{}{}^{12}}+\frac{u^{}}{\xi _0^{}}\mathrm{exp}((r^{}1)/\xi _0^{})$$
(11)
respectively. The same parameters are used as in the monodisperse case.

## III Results

We calculate the radial distribution function to study the stucture of the aggregations. The normalized two dimensional radial distribution function is defined as

$$g_0(r^{})=\frac{1}{8\mathrm{\Phi }N\mathrm{\Delta }}\underset{i}{}\underset{ji}{}\frac{\delta (r^{}r_{ij}^{})}{r_{ij}^{}}$$
(12)
where $`\mathrm{\Delta }`$ is the width of histogram of $`g_0`$. $`g_0`$ is normalised such that it tends to 1 at long distances.

For the analysis of simulations of two types of particles, we introduce partial radial distribution functions, namely, $`g_{AA}`$ is the contribution of $`g_0`$ by radial distances between particles of type A, $`g_{BB}`$ the contributions of type B, and $`g_{AB}`$ the contributions of radial distances between different types A and B. It is clear that $`g_0=g_{AA}+g_{BB}+g_{AB}`$. The distances are counted over the particles in the cells and their periodic images, in the same way as calculating the forces. The calculations are averaged over 20 ensembles for better statistics.

We also calculated the mean coordination number, counted over the particles and their periodic images and averaged over 20 ensembles.

1. Nonequilibrium structure formation for two types of particles, driven by type-dependent forces: We will focus on the time evolution of structure formation as well as on the structural properties of the resulting aggregation by studying the radial distribution function (RDF). We plot the time series of $`g_0`$ of aggregation in Fig.1 with four panels (a)–(d), from initial to final time, both for one type and two types of particles (shown on the same figure). The same figure will be used in describing results in item 2 below.

We plot the partial RDF’s ($`g_{AA}`$, $`g_{BB}`$, $`g_{AB}`$) as well as the total RDF ($`g_0`$) in Fig.2 with six panels (a)–(f), corresponding to the six area fractions studied, for the binary colloidal case. The resulting aggregation is well described by a binary alloy gelation.

2. We will compare the structural properties to those for one type of particles (monodisperse case). We plot the RDF ($`g_0`$) in Fig.3 with six panels (a)–(f), corresponding to the six area fractions studied, for the monodisperse case. From the results, it is evident that there are significant differences between the RDF’s of the two cases. In this regard, we will also study the statistical geometry by computing the average coordination number of the resulting aggregation.

We plot the average coordination number $`Z`$ versus area fraction in Fig.4 both for one type and two types of particles. We find that the average coordination number is generally larger in the monodisperse case for all area fractions considered. It turns out that more close-packed structures (in fact hexagonal) are formed for monodisperse colloidal aggregation.

3. We will also study the partial RDF’s in more detail. By symmetry, the resulting binary colloidal alloy is statistically bipartite, i.e., $`g_{AA}=g_{BB}`$. That means that the resulting lattice can be decomposed into two statistically identical lattice. According to our results in Fig.2, the resulting aggregation can be decomposed into local square order (tetravalent sites) and local hexagonal order (trivalent sites). In the limit of a large attractive force between the different types of particles, however, we may enhance tetravalent sites at the expense of trivalent sites. In this regard, we will do additional simulations for the largest area fraction (corresponding to $`N=350`$), with a larger $`u_1=12`$, other parameters being fixed, so as to suppress the trivalent sites. In Fig.5, we plot the partial RDF’s for the $`u_1=12`$ case. The $`u_1=6`$ case is plotted on the same figure for comparison.

## IV Discussion and conclusion

Here a few comments on our results are in order. The self-assembly of two types of particles can be realized by a recent experimental demonstration of DNA-assisted self-assembly of nanoparticles . While the present investigation has been on isotropic interparticle forces, it is instructive to extend the present work to polydisperse electrorheological (ER) fluids, in which the suspended particles can have different dielectric permittivities . In ER fluids, the polarized particles aggregate under the influence of anisotropic dipolar forces. By tuning the strength of the applied field, it is possible to realize a structure transformation from the body-centered tetragonal to the face-centered cubic structures . Thus, by considering two different types of particles driven both by isotropic and anisotropic interparticle forces, we may obtain a diversity of structures with potential applications in photonic band-gap materials.

In conclusion, we have performed a detailed molecular dynamics simulation to study the nonequilibrium structure formation for two types of particles in a colloidal suspension, driven by type-dependent forces. We examined the time evolution of structure formation as well as on the structural properties of the resulting aggregation by studying the radial distribution function (RDF). The resulting aggregation is well described by a binary colloidal gelation. We find that the average coordination number is generally larger in the monodisperse case for all area fractions considered. Thus, by means of heterogeneous aggregation, it is possible to obtain a wide variety of structures while more close-packed structures are formed for monodisperse colloidal aggregation.

## Acknowledgments

This work was supported by the Research Grants Council of the Hong Kong SAR Government under grant CUHK 4245/01P. K.W.Y. acknowledges the hospitality received during his participation in the Adriatico Research Conference on Interaction and Assembly of Biomolecules, hosted by the International Center for Theoretical Physics at Italy, where the present work was initiated. We acknowledge useful discussion with Dr. Jones T. K. Wan.

Fig.1/Wong and Yu

Fig.2/Wong and Yu

Fig.3/Wong and Yu

Fig.4/Wong and Yu

Fig.5/Wong and Yu


text:
# Maximum Speedup in Quantum Search : 𝑂(1) Running Time

## Acknowledgement

J. Bae is supported in part by the Hanyang University Fellowship and Y. Kwon is supported in part by the Fund of Hanyang University.


text:
# 1 Percolation cumulant of pure-site clusters in the 2D Ising model as a function of 𝛽=𝐽/𝑘𝑇 for different lattice sizes. The curves cross remarkably at the same point, in agreement with the critical point of the magnetization transition (vertical dotted line in the plot).

Site Percolation and Phase Transitions in Two Dimensions

Santo Fortunato

Fakultät für Physik, Universität Bielefeld

D-33501 Bielefeld, Germany

Abstract:

The properties of the pure-site clusters of spin models, i.e. the clusters which are obtained by joining nearest-neighbour spins of the same sign, are here investigated. In the Ising model in two dimensions it is known that such clusters undergo a percolation transition exactly at the critical point. We show that this result is valid for a wide class of bidimensional systems undergoing a continuous magnetization transition. We provide numerical evidence for discrete as well as for continuous spin models, including $`SU(N)`$ lattice gauge theories. The critical percolation exponents do not coincide with the ones of the thermal transition, but they are the same for models belonging to the same universality class.

The idea to explain the mechanism of a phase transition in terms of the interplay of geometrical structures which can be identified in the system is quite old : the growth of such structures close to the transition represents the increase of the range of the correlation between different parts of the system, so that the formation of a spanning structure corresponds to a divergent correlation length and to the transition to a new state of global order for the system.

Here we exclusively refer to lattice models and the geometrical structures of interest are usually clusters of neighbouring particles (spins). Percolation theory  is the ideal framework to deal with clusters. Near the critical point of the percolation transition there is power-law behaviour for the percolation variables and one can define an analogous set of critical exponents as for standard thermal transitions. The main percolation variables are:

* the percolation strength $`P`$, i.e. the probability that a site chosen at random belongs to a percolating cluster;
* the average cluster size $`S`$,

$$S=\frac{_sn_ss^2}{_sn_ss},$$
(1)
where $`n_s`$ is the number of clusters with $`s`$ sites and the sums exclude eventual percolating clusters.

One can in principle define the clusters arbitrarily, but, in order to reproduce the critical behavior of the model, the following conditions must be satisfied:

* the percolation point must coincide with the thermal critical point;
* the connectedness length (average cluster radius) diverges as the thermal correlation length (same exponent);
* the percolation strength $`P`$ near the threshold varies like the order parameter $`m`$ of the model (same exponent);
* the average cluster size $`S`$ diverges as the physical susceptibility $`\chi `$ (same exponent).

The first system to be investigated was of course the Ising model. If one takes an Ising configuration one can immediately isolate the ”classical” magnetic domains, i.e. the clusters formed by binding to each other all nearest-neighbouring sites carrying spins of the same sign. It turns out that these clusters have indeed interesting properties in two dimensions, since their percolation temperature coincides with the critical temperature of the magnetization transition ; the percolation exponents, however, do not coincide with the Ising exponents .

In this paper we show that the result is quite general in two dimensions, being valid for several models undergoing a continuous magnetization transition. In particular we shall see that the presence of antiferromagnetic interactions does not affect the result, at least to the extent of the cases studied here. Moreover the critical percolation exponents do not randomly vary from one model to the other but are uniquely fixed by the universality class the model belongs to.

Our results are based on Monte Carlo simulations on square lattices of various systems near criticality. The models we studied can be divided in two groups: models with $`Z(2)`$ global symmetry and a magnetization transition with Ising exponents and models with $`Z(3)`$ global symmetry and a magnetization transition with exponents belonging to the 2-dimensional 3-state Potts model universality class. The systems belonging to the first group are:

1. the Ising model, $`=J_{ij}s_is_j`$ ($`J>0`$, $`s_i=\pm 1`$);
2. a model with nearest-neighbour (NN) ferromagnetic coupling and a weaker next-to-nearest (NTN) antiferromagnetic coupling: $`=J_1_{NN}s_is_jJ_2_{NTN}s_is_j`$ ($`J_1>0`$, $`J_2<0`$, $`|J_2/J_1|=1/10`$, $`s_i=\pm 1`$);
3. the continuous Ising model, $`=J_{ij}S_iS_j`$ ($`J>0`$, $`1S_i+1`$);
4. SU(2) pure gauge theory in 2+1 dimensions.

The models belonging to the second group are:

1. the 3-state Potts model, $`=J_{ij}\delta (s_i,s_j)`$ ($`J>0`$, $`s_i=1,2,3`$);
2. a model obtained by adding to 1) a weaker next-to-nearest (NTN) antiferromagnetic coupling: $`=J_1_{NN}\delta (s_i,s_j)J_2_{NTN}\delta (s_i,s_j)`$ ($`J_1>0`$, $`J_2<0`$, $`|J_2/J_1|=1/10`$, $`s_i=1,2,3`$);

In each case we produced the thermal equilibrium configurations by using standard Monte Carlo algorithms, like Metropolis or heat bath; whenever we could (Ex. models 1 and 3 of the first group, model 1 of the second) we adopted cluster algorithms in order to reduce the autocorrelation time***For the models 2 of both groups a cluster dynamics exists, but it does not bring great advantages because of the huge size of the clusters to be flipped at the critical point. At each iteration we calculated the energy and the magnetization; the latter is necessary to study the thermal transition. For the identification of the pure-site clusters in each configuration we made use of the algorithm devised by Hoshen and Kopelman , with free boundary conditions. Finally we determined the percolation strength $`P`$, the average cluster size $`S`$ and the size $`S_M`$ of the largest cluster of each configuration, from which we can calculate the fractal dimension $`D`$ of the percolating cluster. We say that a cluster percolates if it spans the lattice from top to bottom. The finite size scaling laws at the critical temperature $`T_p`$ for the percolation variables read

$`P(T_p)`$ $``$ $`L^{\beta _p/\nu _p}`$ (2)
$`S(T_p)`$ $``$ $`L^{\gamma _p/\nu _p}`$ (3)
$`S_M(T_p)`$ $``$ $`L^D,`$ (4)
where $`L`$ is the lattice side and $`\beta _p`$, $`\gamma _p`$, $`\nu _p`$ are critical exponents of the percolation transition. Moreover, to study the percolation transition it is helpful to define also the percolation cumulant. It is the probability of having percolation at a given temperature and lattice size, i.e., the quantity obtained by dividing the number of ”percolating” samples by the total number of analyzed configurations. This variable has three remarkable properties:

* if one plots it as a function of $`T`$, all curves corresponding to different lattice sizes cross at the same temperature $`T_p`$, which marks the threshold of the percolation transition;
* the percolation cumulants for different values of the lattice size $`L`$ coincide, if considered as functions of $`t_pL^{1/\nu _p}`$ ($`t_p=(TT_p)/T_p`$, $`\nu _p`$ is the exponent of the percolation correlation length);
* the value of the percolation cumulant at $`T_p`$ is a universal quantity, i.e. it labels a well defined set of critical indices.

These features allow a rather precise determination of the critical point and a fairly good estimate of the critical exponent $`\nu _p`$.

We start to expose our results from the models with an Ising-like transition. We included the Ising model itself in order to reproduce the results of  and to precisely determine the percolation exponents. We simulated the model on several lattices, from $`600^2`$ to $`2000^2`$. The number of measurements we have taken ranged from a minimum of 20000 for the large $`2000^2`$ lattice to over 100000 for the others. From the crossing of the cumulants (Fig. 1) we obtained for the percolation temperature $`J/kT_p=0.44069(1)`$, in excellent agreement with the critical temperature $`J/kT_c=\mathrm{log}(1+\sqrt{2})/2=0.44068679\mathrm{}`$. For all models we investigated it is possible to perform a precise interpolation of the raw data by applying the density of state method  both to the thermal and to the percolation variables. The curves in Fig. 1 are examples of such interpolations; we did not put the errors to show more clearly the crossing point and the fact that the expected behaviour is already very well represented by the mean values, which is due to the high statistics we could reach in this simple case.

By rescaling the percolation cumulant curves as described above, we obtain Fig. 2 if we choose for the exponent $`\nu _p`$ the value of the thermal 2D Ising exponent $`\nu _{Is}=1`$. This time we have also drawn the errors, and we see that all data relative to different lattice sizes fall onto one and the same curve, which shows that $`\nu _p=\nu _{Is}=1`$, as suggested in .

The percolation exponents and the threshold value of the cumulant are listed in Table I (Model 1), where we also put the critical percolation indices of the other three models of this group. The values of $`\gamma _p`$, $`\beta _p`$ and of the fractal dimension $`D`$ are in good accord with the theoretical predictions of  (first line of Table I).

For the model with competitive interactions we plot the percolation cumulant for two lattice sizes, $`100^2`$ and $`200^2`$ (Fig. 3). Our statistics for this system is of 20000 measurements for each temperature and lattice size, we could not increase it because the algorithms one can use are not very efficient in this case. However, the result is clear: we see that the percolation temperature agrees with the magnetization one (dotted vertical line in the plot). More precisely we find $`J_1/kT_p=0.51422(12)`$, to be compared with $`J_1/kT_c=0.51418(10)`$. The values of the critical percolation exponents agree with the ones of the 2D Ising model as we can see in Table I (Model 2).

The magnetization transition of the continuous Ising model was studied in . The value of the critical temperature is $`J/kT_c=1.09312(1)`$. For our simulations we made use of the same algorithm described in , which consists in a combination of heat bath steps and cluster flippings; this algorithm is quite efficient and we could push our statistics up to 100000 measurements for each lattice size. We remark that in this case we have to do with a continuous spin variable, so that we join nearest-neighbouring like-signed spins, although their absolute values are, in general, different. Fig. 4 shows the percolation cumulant for the four lattice sizes we took: $`100^2`$, $`200^2`$, $`300^2`$, $`400^2`$. Due to our high statistics, we did not draw the errors in our plot, and we see that the curves of the mean values cross remarkably at the same point, which coincides with the thermal critical point (dotted vertical line in the figure). Our estimate of the percolation temperature is $`J/kT_p=1.09311(3)`$.

We notice once again that the critical indices are, within errors, the same we found for the previous two models (Model 3 in Table I).

Finally we considered the more involved case of $`SU(2)`$ lattice gauge theory, in the pure gluonic sector. It is known that in this case there is a confinement-deconfinement transition with Ising exponents. The order parameter is the lattice average $`L`$ of the Polyakov loop, which is a continuous variable like the spins of the continuous Ising model we examined above. We analyzed the 2+1 dimensional case, which corresponds to two space dimensions. The number of lattice spacings in the imaginary time direction is $`N_\tau =2`$. We took four lattices: $`64^2`$, $`96^2`$, $`128^2`$ and $`200^2`$. The statistics ranged from 20000 to 50000 measurements. Our value for the critical coupling is $`\beta _c=4/g_{c}^{}{}_{}{}^{2}=3.4504(11)`$, which is by an order of magnitude more precise than the value reported in . The original $`SU(2)`$ matrix configurations on the links of the 2+1 dimensional lattice are projected onto Polyakov loop configurations on the sites of a 2-dimensional lattice. For the percolation analysis we investigated the Polyakov loop configurations. In order to build the clusters we proceeded as for the continuous Ising model, by joining nearest-neighbouring sites carrying like-signed values of the Polyakov loop variable. In Fig. 5 we plotted the percolation cumulant curves as a function of the coupling $`\beta `$ for two lattice sizes, $`96^2`$ and $`200^2`$. Also here the coincidence of the percolation with the critical point seems to be clear. The percolation coupling we determined is $`\beta _p=4/g_{p}^{}{}_{}{}^{2}=3.4501(14)`$. The critical percolation indices are again the same we found so far (Model 4 in Table I).

We complete our investigation by studying the models of the second group. Now we have three spin states and the clusters are simply obtained by binding nearest-neighbouring sites in the same spin state. The fact that the pure-site clusters percolate at $`T_c`$ for the 3-state Potts model in two dimensions was found, though not rigorously, in . The percolation temperature that we determined is $`J/kT_p=1.00511(9)`$, in agreement with the thermal threshold $`J/kT_c=\mathrm{log}(1+\sqrt{3})=1.00505\mathrm{}`$. The exponents of the percolation transition are listed in Table II (Model 1). We see that they are different from the thermal exponents ($`\beta /\nu =2/15`$, $`\gamma /\nu =26/15`$, $`\nu =5/6`$), except $`\nu _p`$, which, like in the Ising case, coincides with the thermal value $`5/6`$. We remark that the values of the critical indices are in very good accord with the predictions of  (first line of Table II)). The percolation exponents are also different from the ones we have found for the models in the 2D Ising universality class.

The analysis of the model with competitive interactions leads to the same conclusions: the percolation point is $`J_1/kT_p=1.1670(7)`$, in accord with the magnetization temperature $`J_1/kT_c=1.1665(9)`$. In Fig. 6 we plotted both the Binder and the percolation cumulant, so that we can have a visual comparison of the two thresholds. The critical percolation indices agree with the ones we have found for the 3-state Potts model, as we can see in Table II (Model 2).

In conclusion, we have found that the percolation point of the pure-site clusters of various models coincides with the critical point of the thermal magnetization transition. We examined theories with discrete and continuous spins, involving interactions beyond the fundamental nearest-neighbour one, and we included competitive interactions by means of antiferromagnetic couplings. In particular, we stress that the effective theory of the Polyakov loop for $`SU(2)`$ gauge theory (and in general for $`SU(N)`$) includes many couplings, not only nearest-neighbour spin-spin interactionsThe couplings can involve spins which are far from each other, a higher number of spins (plaquette interactions, six-spins couplings, etc.), and also self-interactions., and that such couplings are ferromagnetic but also antiferromagnetic (see ). Therefore, the fact that our result is valid even for $`SU(2)`$ suggests that it is a rather general feature of bidimensional models with a continuous magnetization transition. It would be interesting to check whether this property of the pure-site clusters always holds in two dimensions. We remark as well that, even if we analyzed models with magnetization transitions, the result is probably also true for systems undergoing a phase transition towards antiferromagnetic ground states, since in this case the definition of the clusters can be trivially extended by joining nearest-neighbouring spins of opposite signsIn the Ising model the equivalence of the two cases is evident..

Moreover, we found that the critical percolation indices are the same for systems in the same universality class, although they do not coincide with the thermal indices. This shows that, at criticality, the pure-site clusters of all models in the same universality class are, virtually, indistinguishable from each other as far as their size distributions are concerned (although they could differ in their topological properties). So, in many respects, the behaviour of the pure-site clusters of a model is as universal as other features of the model, e.g. its critical exponents. That indicates that there is a close relationship between such clusters and the critical behaviour of the system. It is indeed possible to show that the pure-site clusters of a model ”contain” the clusters which reproduce its critical behaviour in the sense explained at the beginning of this paper .

The fact that the pure-site clusters may percolate at the critical point also for systems different from the Ising model was suggested in ; however, this prediction concerned just models with ferromagnetic interactions and is a conjecture. This work extends the result to models with competitive interactions to which the theorems and arguments used in  cannot be applied. We guess that the property we have illustrated here is due to the fact that the pure-site clusters of the Ising model, for instance, are quite compact objects. This is confirmed  by the fact that, in two dimensions, even if one randomly breaks the bonds between nearest-neighbour spins of the same sign with a probability $`p_B=\mathrm{exp}(2J/kT)`$, the corresponding site-bond clusters keep percolating at the critical temperature $`T_c`$. That means that the infinite pure-site cluster at $`T_c`$ remains an infinite cluster even if we randomly break many of the bonds between the spins belonging to it, which considerably reduces its size. The presence of antiferromagnetic couplings also acts like a bond-breaking factor, which tends to reduce the size of the original pure-site clusters because of the fact that they favour the presence of anti-aligned spins in the system. Nevertheless, if the anti-alignement interactions are weaker that the main nearest-neighbour ferromagnetic coupling, it is likely that they do not succeed in breaking the original spanning structure into finite clusters. That could also explain the fact that the critical percolation exponents are the same as in the original model, since we can imagine that the antiferromagnetic interactions do not perturb enough the renormalization group trajectories in the (percolation) coupling space, so that they would keep on converging to the same fixed point, exactly as it happens for the thermal transition. That naturally leads to the question: how small should the antiferromagnetic couplings be compared to the fundamental ferromagnetic one, so that the pure-site clusters still percolate at $`T_c`$? For the systems we have studied the ratio was of 1:10; we have evidence that the result holds also for a ratio 3:10, we did not present it here because the statistics is quite low. It would be particularly interesting to check whether the property holds as long as the competition between the different interactions still allows a phase transition with spin-ordering to take place.

I would like to thank my PhD advisor, professor Helmut Satz, for introducing me in this exciting field and for his support over the years. I also thank professor A. Coniglio for helpful discussions. I gratefully acknowledge the financial support of the TMR network ERBFMRX-CT-970122 and the DFG Forschergruppe FOR 339/1-2.