text:
# Collective Origin of the Coexistence of Apparent RMT Noise and Factors in Large Sample Correlation Matrices

## Abstract

Through simple analytical calculations and numerical simulations, we demonstrate the generic existence of a self-organized macroscopic state in any large multivariate system possessing non-vanishing average correlations between a finite fraction of all pairs of elements. The coexistence of an eigenvalue spectrum predicted by random matrix theory (RMT) and a few very large eigenvalues in large empirical correlation matrices is shown to result from a bottom-up collective effect of the underlying time series rather than a top-down impact of factors. Our results, in excellent agreement with previous results obtained on large financial correlation matrices, show that there is relevant information also in the bulk of the eigenvalue spectrum and rationalize the presence of market factors previously introduced in an ad hoc manner.

Since Wigner’s seminal idea to apply random matrix theory (RMT) to interpret the complex spectrum of energy levels in nuclear physics Wigner , RMT has made enormous progress Mehta  with many applications in physical sciences and elsewhere such as in meteorology SP01  and image processing SM . A new application was proposed a few years ago to the problem of correlations between financial assets and to the portfolio optimization problem. It was shown that, among the eigenvalues and principal components of the empirical correlation matrix of the returns of hundreds of asset on the New York Stock Exchange (NYSE), apart from the few highest eigenvalues, the marginal distribution of the other eigenvalues and eigenvectors closely resembles the spectral distribution of a positive symmetric random matrix with maximum entropy, suggesting that the correlation matrix does not contain any specific information beyond these few largest eigenvalues and eigenvectors Laloux . These results apparently invalidate the standard mean-variance portfolio optimization theory Marco  consecrated by the financial industry riskmetrics  and seemingly support the rationale behind factor models such as the capital asset pricing model (CAPM) CAPM  and the arbitrage pricing theory (APT) APT , where the correlations between a large number of assets are represented through a small number of so-called market factors. Indeed, if the spectrum of eigenvalues of the empirical covariance or correlation matrices are predicted by RMT, it seems natural to conclude that there is no usable information in these matrices and that empirical covariance matrices should not be used for portfolio optimization. In contrast, if one detects deviations between the universal – and therefore non-informative – part of the spectral properties of empirically estimated covariance and correlation matrices and those of the relevant ensemble of random matrices, this may quantify the amount of real information that can be used in portfolio optimization from the “noise” that should be discarded.

More generally, in many different scientific fields, one needs to determine the nature and amount of information contained in large covariance and correlation matrices. This occurs as soon as one attempts to estimate very large covariance and correlation matrices in multivariate dynamics of systems exhibiting non-Gaussian fluctuations with fat tails and/or long-range time correlations with intermittency. In such cases, the convergence of the estimators of the large covariance and correlation matrices is often too slow for all practical purposes. The problem becomes even more complex with time-varying variances and covariances as occurs in systems with heteroskedasticity Engle  or with regime-switching regimeswitch . A prominent example where such difficulties arise is the data-assimilation problem in engineering and in meteorology where forecasting is combined with observations iteratively through the Kalman filter, based on the estimation and forward prediction of large covariance matrices Kalman .

As we said in the context of financial time series, the rescuing strategy is to invoke the existence of a few dominant factors, such as an overall market factor and the factors related to firm size, firm industry and book-to-market equity, thought to embody most of the relevant dependence structure between the studied time series fama . Indeed, there is no doubt that observed equity prices respond to a wide variety of unanticipated factors, but there is much weaker evidence that expected returns are higher for equities that are more sensitive to these factors, as required by Markowitz’s mean-variance theory, by the CAPM and the APT Roll . This severe failure of the most fundamental finance theories could conceivably be attributable to an inappropriate proxy for the market portfolio, but nobody has been able to show that this is really the correct explanation. This remark constitutes the crux of the problem: the factors invoked to model the cross-sectional dependence between assets are not known in general and are either postulated based on economic intuition in financial studies or obtained as black box results in the recent analyses using RMT Laloux .

Here, we show that the existence of factors results from a collective effect of the assets, similar to the emergence of a macroscopic self-organization of interacting microscopic constituents. For this, we unravel the general physical origin of the large eigenvalues of large covariance and correlation matrices and provide a complete understanding of the coexistence of features resembling properties of random matrices and of large “anomalous” eigenvalues. Through simple analytical calculations and numerical simulations, we demonstrate the generic existence of a self-organized macroscopic state in any large system possessing non-vanishing average correlations between a finite fraction of all pairs of elements.

Let us first consider a large system of size $`N`$ with correlation matrix $`C`$ in which every non-diagonal pairs of elements exhibits the same correlation coefficient $`C_{ij}=\rho `$ for $`ij`$ and $`C_{ii}=1`$. Its eigenvalues are

$$\lambda _1=1+(N1)\rho \mathrm{and}\lambda _{i2}=1\rho $$
(1)
with multiplicity $`N1`$ and with $`\rho (0,1)`$ in order for the correlation matrix to remain positive definite. Thus, in the thermodynamics limit $`N\mathrm{}`$, even for a weak positive correlation $`\rho 0`$ (with $`\rho N1`$), a very large eigenvalue appears, associated with the delocalized eigenvector $`v_1=(1/\sqrt{N})(1,1,\mathrm{},1)`$, which dominates completely the correlation structure of the system. This trivial example stresses that the key point for the emergence of a large eigenvalue is not the strength of the correlations, provided that they do not vanish, but the large size $`N`$ of the system.

This result (1) still holds qualitatively when the correlation coefficients are all distinct. To see this, it is convenient to use a perturbation approach. We thus add a small random component to each correlation coefficient:

$$C_{ij}=\rho +ϵa_{ij}\mathrm{for}ij,$$
(2)
where the coefficients $`a_{ij}=a_{ji}`$ have zero mean, variance $`\sigma ^2`$ and are independently distributed (There are additional constraints on the support of the distribution of the $`a_{ij}`$’s in order for the matrix $`C_{ij}`$ to remain positive definite with probability one). The determination of the eigenvalues and eigenfunctions of $`C_{ij}`$ is performed using the perturbation theory developed in quantum mechanics cohen  up to the second order in $`ϵ`$. We find that the largest eigenvalue becomes

$$\mathrm{E}[\lambda _1]=(N1)\rho +1+\frac{(N1)(N2)}{N^2}\frac{ϵ^2\sigma ^2}{\rho }+𝒪(ϵ^3)$$
(3)
while, at the same order, the corresponding eigenvector $`v_1`$ remains unchanged. The degeneracy of the eigenvalue $`\lambda =1\rho `$ is broken and leads to a complex set of smaller eigenvalues described below.

In fact, this result (3) can be generalized to the non-perturbative domain of any correlation matrix with independent random coefficients $`C_{ij}`$, provided that they have the same mean value $`\rho `$ and variance $`\sigma ^2`$. Indeed, it has been shown FK81  that, in such a case, the expectations of the largest and second largest eigenvalues are

$`\mathrm{E}[\lambda _1]`$ $`=`$ $`(N1)\rho +1+\sigma ^2/\rho +o(1),`$ (4)
$`\mathrm{E}[\lambda _2]`$ $``$ $`2\sigma \sqrt{N}+𝒪(N^{1/3}\mathrm{log}N).`$ (5)
Moreover, the statistical fluctuations of these two largest eigenvalues are asymptotically (for large fluctuations $`t>𝒪(\sqrt{N})`$) bounded by a Gaussian distribution according to the following large deviation theorem

$$\mathrm{Pr}\{|\lambda _{1,2}\mathrm{E}[\lambda _{1,2}]|t\}e^{c_{1,2}t^2},$$
(6)
for some positive constant $`c_{1,2}`$ KV00 .

This result is very different from that obtained when the mean value $`\rho `$ vanishes. In such a case, the distribution of eigenvalues of the random matrix $`C`$ is given by the semi-circle law Mehta . However, due to the presence of the ones on the main diagonal of the correlation matrix $`C`$, the center of the circle is not at the origin but at the point $`\lambda =1`$. Thus, the distribution of the eigenvalues of random correlation matrices with zero mean correlation coefficients is a semi-circle of radius $`2\sigma \sqrt{N}`$ centered at $`\lambda =1`$.

The result (4) is deeply related to the so-called “friendship theorem” in mathematical graph theory, which states that, in any finite graph such that any two vertices have exactly one common neighbor, there is one and only one vertex adjacent to all other vertices Erdos . A more heuristic but equivalent statement is that, in a group of people such that any pair of persons have exactly one common friend, there is always one person (the “politician”) who is the friend of everybody. The connection is established by taking the non-diagonal entries $`C_{ij}`$ ($`ij`$) equal to Bernouilli random variable with parameter $`\rho `$, that is, $`Pr[C_{ij}=1]=\rho `$ and $`Pr[C_{ij}=0]=1\rho `$. Then, the matrix $`C_{ij}I`$, where $`I`$ is the unit matrix, becomes nothing but the adjacency matrix of the random graph $`G(N,\rho )`$ KV00 . The proof of Erdos  of the “friendship theorem” indeed relies on the $`N`$-dependence of the largest eigenvalue and on the $`\sqrt{N}`$-dependence of the second largest eigenvalue of $`C_{ij}`$ as given by (4) and (5).

Figure 1 shows the distribution of eigenvalues of a random correlation matrix. The inset shows the largest eigenvalue lying at the predicting size $`\rho N=56.8`$, while the bulk of the eigenvalues are much smaller and are described by a modified semi-circle law centered on $`\lambda =1\rho `$, in the limit of large $`N`$. The result on the largest eigenvalue emerging from the collective effect of the cross-correlation between all $`N(N1)/2`$ pairs provides a novel perspective to the observation Rollcrash  that the only reasonable explanation for the simultaneous crash of 23 stock markets worldwide in October 1987 is the impact of a world market factor: according to our demonstration, the simultaneous occurrence of significant correlations between the markets worldwide is bound to lead to the existence of an extremely large eigenvalue, the world market factor constructed by … a linear combination of the 23 stock markets! What our result shows is that invoking factors to explain the cross-sectional structure of stock returns is cursed by the chicken-and-egg problem: factors exist because stocks are correlated; stocks are correlated because of common factors impacting them.

Figure 2 shows the eigenvalues distribution of the sample correlation matrix reconstructed by sampling $`N=406`$ time series of length $`T=1309`$ generated with a given correlation matrix $`C`$ with theoretical spectrum shown in figure 1. The largest eigenvalue is again very close to the prediction $`\rho N=56.8`$ while the bulk of the distribution departs very strongly from the semi-circle law and is not far from the Wishart prediction, as expected from the definition of the Wishart ensemble as the ensemble of sample covariance matrices of Gaussian distributed time series with unit variance and zero mean. A Kolmogorov test shows however that the bulk of the spectrum is not in the Wishart class, in contradiction with previous claims lacking formal statistical tests Laloux . This result holds for different simulations of the sample correlation matrix and different realizations of the theoretical correlation matrix with the same parameters ($`\rho `$, $`\sigma `$). The statistically significant departure from the Wishart prediction implies that there is actually some information in the bulk of the spectrum of eigenvalues, which is intimately coupled with the existence of the largest eigenvalue. We have also checked that these results remain robust for non-Gaussian distribution of returns as long as the second moments exist. Indeed, correlated time series with multivariate Gaussian or Student distributions with three degrees of freedom (which provide more acceptable proxies for financial time series Gopikrishnan ) give no discernible differences in the spectrum of eigenvalues. This is surprising as the estimator of a correlation coefficient is asymptotically Gaussian for time series with finite fourth moment and Lévy stable otherwise MS2 .

Up to now, we have focused on the collective mechanism at the origin of the very large eigenvalue of order $`N`$. Empirically Laloux , a few other eigenvalues have an amplitude of the order of $`510`$ that deviate significantly from the bulk of the distribution. These eigenvalues cannot be obtained by a matrix of the form (2) with identically independently distributed coefficients $`a_{ij}`$’s. Our analysis provides a very simple constructive mechanism for them. The solution consists in considering, as a first approximation, the block diagonal matrix $`C^{}`$ with diagonal elements made of the matrices $`A_1,\mathrm{},A_p`$ of sizes $`N_1,\mathrm{},N_p`$ with $`N_i=N`$, constructed according to (2) such that each matrix $`A_i`$ has the average correlation coefficient $`\rho _i`$. When the coefficients of the matrix $`C^{}`$ outside the matrices $`A_i`$ are zero, the spectrum of $`C^{}`$ is given by the union of all the spectra of the $`A_i`$’s, which are each dominated by a large eigenvalue $`\lambda _{1,i}\rho _iN_i`$. The spectrum of $`C^{}`$ then exhibits $`p`$ large eigenvalues.

Each block $`A_i`$ can be interpreted as a sector of the economy, including all the companies belonging to a same industrial branch and the eigenvector associated with each largest eigenvalue represents the main factor driving this sector of activity Mantegnapath . For similar sector sizes $`N_i`$ and average correlation coefficients $`\rho _i`$, the largest eigenvalues are of the same order of magnitude. In order to recover a very large unique eigenvalue, we reintroduce some coupling constants outside the block diagonal matrices. A well-known result of perturbation theory in quantum mechanics states that such coupling leads to a repulsion between the eigenstates, which can be observed in figure 3 where $`C^{}`$ has been constructed with three block matrices $`A_1`$, $`A_2`$ and $`A_3`$ and non-zero off-diagonal coupling described in the figure caption. These values allow us to quantitatively replicate the empirical finding of Laloux et al. in Laloux , where the three first eigenvalues are approximately $`\lambda _157`$, $`\lambda _210`$ and $`\lambda _38`$. The bulk of the spectrum (which excludes the three largest eigenvalues) is similar to the Wishart distribution but again statistically different from it as tested with a Kolmogorov test. There is thus significant deviation from the predictions of RMT not only for the largest eigenvalues but also in the bulk.

As a final remark, expressions (3,4) and our numerical tests for a large variety of correlation matrices show that the delocalized eigenvector $`v_1=(1/\sqrt{N})(1,1,\mathrm{},1)`$, associated with the largest eigenvalue is extremely robust and remains (on average) the same for any large system. Thus, even for time-varying correlation matrices – as in finance with important heteroskedastic effects – the composition of the main factor remains almost the same. This can be seen as a generalized limit theorem reflecting the bottom-up organization of broadly correlated time series.

This work was partially supported by the James S. Mc Donnell Foundation 21st century scientist award/studying complex system.


format: text

text:
# “Hidden” Seyfert 2 Galaxies and the X-ray Background

## 1 Introduction

Identification of the classes of objects responsible for the cosmic X-ray background radiation (XRB) has been one of the fundamental pursuits of X-ray astronomy for the past several decades. It now appears that different populations produce the majority of the background in different X-ray energy bands. In the soft 0.5–2 keV range, active galactic nuclei (AGNs) with broad optical emission lines (i.e., Seyfert 1 galaxies and quasars) are the dominant contributors to the XRB (Schmidt et al. 1998). These objects, however, have steep X-ray spectra and cannot explain the flat slope of the XRB spectrum at higher energies. Type 2 (narrow-line) Seyfert galaxies have thus emerged as the most promising candidates for the origin of the hard (2–10 keV) XRB. The soft X-ray fluxes of these objects are heavily absorbed by dense circumnuclear gas, presumably the same material that obscures their broad emission-line regions. Models (Madau, Ghisellini, & Fabian 1994; Comastri et al. 1995; Gilli, Salvati, & Hasinger 2001) and X-ray observations of nearby sources (Moran et al. 2001) have demonstrated that Seyfert 2 galaxies have both the space density and X-ray spectral properties necessary to account for the intensity and spectrum of the hard XRB.

The Chandra X-ray Observatory, which resolves the majority of the XRB in deep exposures, provides an opportunity to confirm the Seyfert 2 hypothesis directly. Surprisingly, however, follow-up optical observations of faint sources in the deepest Chandra images have revealed a significant population of apparently normal galaxies and far fewer Seyfert 2s than expected (e.g., Mushotzky et al. 2000; Barger et al. 2001a, 2001b, 2002; Hornschemeier et al. 2001). “Normal” galaxies represent $``$ 40–60% of the hard X-ray sources classified in deep surveys; many exhibit evidence of absorption in the X-ray band, but emission lines are very weak in their starlight-dominated optical spectra. The multiwavelength properties of the prototypical example of this class are described by Comastri et al. (2002a). The prevailing interpretation is that optically normal, X-ray–luminous galaxies are common at modest redshifts, perhaps reflecting a tendency for AGNs in the past to be more absorbed or have weaker ionizing continua than present-day Seyfert galaxies (Barger et al. 2001a). This hints at the possibility that the X-ray galaxy population has evolved with cosmic time.

Alternatively, we propose that the results of the deep Chandra surveys may be heavily influenced by the limitations of ground-based observing. In spectroscopic observations of nearby AGNs, light is collected through a small aperture centered on the nucleus, which excludes most of the starlight from the host galaxy. The faint sources that produce the hard XRB, in contrast, are very distant (the majority are at redshifts $`z>0.5`$); thus, when they are observed, the entire galaxy (or a large fraction of it) falls within the spectrograph slit. (At $`z1`$, 10 kpc corresponds to $``$ 1$`\stackrel{}{\mathrm{.}}`$5.) Combined with the low spectral resolution that has been employed to date ($``$ 12–20 Å) and the modest signal-to-noise ratios frequently obtained (the optical counterparts are generally quite faint), the additional galaxy light (from stars and H II regions) could lead to the appearance that some distant Chandra sources are associated with normal galaxies rather than with Seyfert 2s.

To investigate this possibility we have obtained integrated spectra of nearby Seyfert 2 galaxies that are known to be absorbed X-ray sources. Our techniques simulate ground-based spectroscopic observations of distant X-ray galaxies, allowing us to evaluate whether the emission-line signatures of their activity can be overwhelmed (or “hidden”) in spectra of their integrated light. This project draws inspiration from a series of papers by Kennicutt (1992a, 1992b), who obtained integrated spectra of nearby galaxies of all types to assist “the classification of unresolved galaxies at large distances \[and\] the spectroscopic identification of extended sources detected in infrared, radio, and X-ray surveys” (Kennicutt 1992a). Our similarly motivated work focuses on the class of objects believed to be responsible for the hard X-ray background.

## 2 Observations and Analysis

The objects for this study are drawn from the distance-limited sample of Seyfert galaxies compiled by Ulvestad & Wilson (1989). There are 31 confirmed Seyfert 2s in the sample (Moran et al. 2000), almost all of which have been observed in the 1–10 keV band with the ASCA satellite (Moran et al. 2001). The Ulvestad & Wilson sample is ideal for this project for several reasons. First, because of its distance-limited nature, the sample is relatively unbiased. Second, the hard X-ray luminosities of the objects ($``$ $`10^{41}`$$`10^{43}`$ ergs s-1) are comparable to those of the “normal” galaxies in deep Chandra surveys (Barger et al. 2001a, 2001b; Hornschemeier et al. 2001). And finally, because all of the objects are located within $``$ 60 Mpc, most are close enough that we can isolate the emission from their nuclei yet distant enough (i.e., with sufficiently small angular diameters) that we can measure their integrated spectra with relative ease.

Our data were acquired with the 3-m Shane reflector at Lick Observatory over the course of three runs in August and November 2001 and February 2002. We employed 600 line mm-1 dispersing elements on each arm of the Kast double spectrograph (Miller & Stone 1993), which provided coverage of the 3400–8100 Å range. For each galaxy, we obtained a 300 s exposure of the nucleus using a $`2^{\prime \prime }`$ slit oriented at the parallactic angle (Filippenko 1982). The nuclear spectra represent $`4^{\prime \prime }`$$`5^{\prime \prime }`$ extractions along the slit and have a resolution of 5–6 Å (FWHM). Integrated spectra of each object were obtained by scanning an east-west slit across the galaxy in declination at a constant rate. Six scans were obtained for most galaxies, although a few of the brighter objects were scanned three to five times. For the scans, the width of the slit was increased to $`6^{\prime \prime }`$ to improve our observing efficiency and to degrade the spectral resolution to 15–16 Å, comparable to that used in observations of faint Chandra sources (Barger et al. 2002; Hornschemeier et al. 2001). The angular distances of the scans ($`30^{\prime \prime }`$$`120^{\prime \prime }`$ in declination) were dictated by the diameters of the galaxies in high-contrast images from the Digitized Sky Survey (DSS). Each scan had a duration of 600 s; combined with the slit width and scan distances, this yielded effective exposures of 30 s to 120 s per scan for each object’s integrated spectrum. These spectra include most of the light from each galaxy above the sky background (i.e., $`35^{\prime \prime }`$ to $`90^{\prime \prime }`$ extractions).

We observed a total of 18 objects from our sample — all those at $`\delta >23^{}`$ with projected angular sizes in right ascension (again, estimated from DSS images) of less than $`2^{}`$. The latter restriction ensured that we would be able to perform accurate sky background subtraction with light that entered the slit near the two ends. The galaxies observed and some pertinent data are listed in Table 1. The optical and X-ray properties of the 13 objects we did not observe do not differ significantly.

## 3 Integrated Spectra of Seyfert 2 Galaxies

Recently published spectroscopic studies of faint Chandra sources are based mainly on data acquired with 8–10 m telescopes (Barger et al. 2002; Hornschemeier et al. 2001; Tozzi et al. 2001). In most respects, the integrated spectra we have obtained accurately simulate observations of distant type 2 AGNs with these instruments. The scan distances and extractions employed for our sample include most of the light from the galaxies; they correspond to physical sizes of $``$ 8–20 kpc, similar to that which is typically included within a 1.5–$`2^{\prime \prime }`$ slit for a distant Chandra source. Our spectral resolution, as discussed above, is comparable to that used in Chandra source surveys. In addition, the effective exposure times associated with our scans are equivalent to 1.5–2 hr Keck exposures of objects $``$ 7–8 mag fainter, i.e., objects with visual magnitudes in the $`m`$ 20–22 range. These roughly correspond to the Keck exposure times and magnitudes of distant Chandra sources that have been classified spectroscopically (Barger et al. 2001a, 2001b, 2002).

Examples of our results are displayed in Figure 1. The nuclear spectrum of NGC 262 (= Mrk 348; lower trace of Fig. 1a) bears all the characteristics of a type 2 Seyfert galaxy: strong \[O II\] $`\lambda 3727`$, \[Ne III\] $`\lambda 3869`$, H$`\beta `$ $`\lambda 4861`$, \[O III\] $`\lambda \lambda 4959,5007`$, \[O I\] $`\lambda 6300`$, H$`\alpha `$ $`\lambda 6563`$, \[N II\] $`\lambda 6584`$, and \[S II\] $`\lambda \lambda 6717,6731`$ emission lines with high \[O III\]/H$`\beta `$, \[O I\]/H$`\alpha `$, \[N II\]/H$`\alpha `$, and \[S II\]/H$`\alpha `$ flux ratios (Veilleux & Osterbrock 1987). These same features are clearly present in the integrated spectrum of the galaxy (upper trace of Fig. 1a), despite the presence of a significant amount of additional starlight. Thus, NGC 262 retains its Seyfert 2 classification in the integrated spectrum.

The spectra in Figures 1b–e tell a much different story. Once again, the nuclear spectra of these galaxies exhibit Seyfert 2 emission-line signatures. In the integrated spectra, however, the lines are far less prominent, and some are completely obliterated. The strongest spectral features are the Ca II $`\lambda \lambda 3934,3968`$, G band ($`\lambda 4300`$), Mg I $`b`$ $`\lambda 5176`$, and Na I D $`\lambda 5893`$ stellar absorption lines typically observed in the spectra of normal early-type galaxies. Note that the \[Ne III\] $`\lambda 3869`$ emission line, which has been used in Chandra source surveys to identify AGNs (Barger et al. 2001a, 2001b; Hornschemeier et al. 2001), is not visible in the integrated spectra. Thus, \[Ne III\] emission cannot be relied upon to indicate the presence of an active nucleus.

The data for NGC 3982 (Fig. 1f) provide an interesting third alternative. Although the emission lines in both the nuclear and integrated spectra of this object are strong, the line-flux ratios in the two differ dramatically. For example, in the nuclear spectrum, \[N II\] $`\lambda 6584`$/H$`\alpha `$ $``$ 1 and \[O III\] $`\lambda 5007`$/H$`\beta `$ $``$ 16, which is typical for a Seyfert 2. In the integrated spectrum, however, \[N II\]/H$`\alpha `$ = 0.37 and \[O III\]/H$`\beta `$ = 0.82, consistent with the values observed in H II regions (e.g., Veilleux & Osterbrock 1987). The \[O I\]/H$`\alpha `$ and \[S II\]/H$`\alpha `$ ratios also resemble those of H II regions. Given only the integrated spectrum of NGC 3982, we would classify it as a starburst galaxy rather than a type 2 Seyfert! We obtain a similar result for another object in our sample, NGC 1667. Thus, it is possible for the emission from star-forming regions to drown out an active nucleus in a galaxy, which sheds light on the nature of the “composite” starburst/Seyfert galaxies uncovered during optical surveys of ROSAT sources (Moran, Halpern, & Helfand 1996).

## 4 Discussion and Conclusions

The last column of Table 1 indicates our assessment of the classifications of our galaxies based on their integrated spectra. While seven objects with strong lines in the nucleus remain as Seyfert 2s in the integrated spectra, 11 others (or $``$ 60%) would not be classified as AGNs. This is comparable to the fraction of Chandra sources ($``$ 40–60%) identified with “normal” galaxies in deep surveys (Barger et al. 2001a, 2001b, 2002; Hornschemeier et al. 2001). Among the normal-looking subset of our sample, nine objects resemble early-type galaxies with red continua and strong absorption features, and two objects have spectra reminiscent of starburst (or post-starburst) galaxies.

All of the objects we observed display some evidence of emission lines in their integrated spectra. While the same is true of all of the “normal” galaxies whose spectra are presented by Barger et al. (2001a, 2001b) and Hornschemeier et al. (2001), 19% of the sources classified by Barger et al. (2002) are described as “absorption-line” galaxies. Some of these objects may in fact have faint emission lines — weak \[O II\] and \[O III\] emission is visible in the average absorption-line galaxy spectrum shown in Figure 5h of Barger et al. (2002). Others, however, may be devoid of emission features in the data obtained. Two factors are likely to contribute to the relatively high incidence of absorption-line systems in the Chandra surveys:

(1) The H$`\alpha `$ and \[N II\] emission lines can be the strongest lines in the spectra of Seyfert galaxies; as such, they are vital for accurate spectroscopic classifications. In fact, the classifications of two of the persistent Seyfert 2s in our sample (Mrk 1066 and MCG $`+`$01-27-020) would be ambiguous without information about the lines near H$`\alpha `$. However, $``$ 70% of the spectra of supposedly normal galaxies published by Barger et al. and Hornschemeier et al. do not cover the H$`\alpha `$ region because of the high redshifts of the objects. It is possible that some of the absorption-line galaxies would exhibit emission lines if this region were observed.

(2) Our relatively small sample of nearby objects does not include all possible types of host galaxies. In particular, galaxies with high optical luminosities are absent. Using the expression given in § 12.5 of Barger et al. (2002) and the information listed in Table 1, we find that 14 (78%) of the objects we observed have $`B`$-band luminosities below $`10^{43}`$ ergs s-1, and all have $`L_B<2\times 10^{43}`$ ergs s-1. In contrast, the majority of the sources spectroscopically classified by Barger et al. have $`L_B`$ in excess of $`10^{43}`$ ergs s-1. Compared to the objects in our sample, the integrated light of these high-$`L_B`$ galaxies could more easily overwhelm the emission from moderately luminous Seyfert nulcei. Host-galaxy luminosity may also be related to the detection of optically dull, X-ray–bright objects at fairly low redshifts ($`z<0.4`$; see Comastri et al. 2002b). The possibility that some host-galaxy light is excluded in the spectra of these objects is probably offset by the fact that their host galaxies are quite luminous — nine of the ten “normal” galaxies in the Comastri et al. (2002b) study have $`L_B`$ greater than a few times $`10^{43}`$ ergs s-1.

Thus, as we have speculated previously (Moran et al. 2001), it is possible to “hide” the true nature of classical type 2 AGNs in their integrated optical spectra. Low signal-to-noise ratio and inadequate spectral coverage due to redshift effects can hinder the classification process further. We conclude, therefore, that the demographics of the distant X-ray galaxy population under investigation in deep Chandra surveys are more likely to be influenced by the limitations of ground-based spectroscopic observations than by some evolution of the population. The majority of the apparently normal galaxies that have turned up in these surveys, at least those with flat (presumably absorbed) X-ray spectra and/or substantial X-ray luminosities, would probably have optical spectra similar to those of nearby Seyfert 2s if their nuclei could be isolated spatially. If so, this would verify that Seyfert 2 galaxies are the dominant contributors to the hard XRB, as prior research has suggested. Host galaxy dilution is likely to impact any ground-based spectroscopic survey of distant galaxies, leading to underestimates of the numbers of narrow-line AGNs in such surveys.

We would like to express our thanks to Bill Holzapfel for generously donating the August 2001 Lick time that enabled us to get this project underway. This work was supported by NASA grant NAG5-3556 to A. V. F.


text:
# Managing catastrophic changes in a collective

## 1 Introduction

Understanding the relationship between the overall macroscopic performance of a collective, and its design at the microscopic level, is a high priority for the collectives field (see, for example, Ref. ). Typically such performance can be measured in terms of a macroscopic variable of the system which fluctuates in time: for example the level of wastage of a global resource. An example which can be seen in other Chapters of this book, is the discussion of the fluctuations in attendance in the El Farol bar game and related multi-agent models such as the Minority Game. In those cases, one is particularly interested in understanding how to minimize the typical size of such fluctuations and hence ‘optimize’ performance by minimizing the average wastage. In practice, however, typical optimization schemes which focus on minimizing the variance of some global quantity may not be the most relevant. Biology is a wonderful example of a non-optimal system, yet one which does a fantastic job of avoiding catastrophic large changes. The human body, for example, has a complex network of feedback loops set in place in order to ‘survive’. While very few of us exist at our peak level, and even fewer of us hold world records in athletics, we do a pretty good job of coping with everything our changing bodies and changing environment can throw at us. Moreover, the human system can self-manage this difficult situation for up to 100 years, an unimaginable feat for any computer-based system. In short, we may not be optimal in the sense of minimizing wastage, yet our system can handle sudden changes in the environment (i.e. exogenous changes) and also typically manages to resist the tendency to self-generate large unexpected (endogenous) changes or ‘system crashes’.

Large unexpected changes, or so-called extreme events, happen infrequently, yet tend to dictate the long-term dynamical behaviour of real-world systems in disciplines as diverse as biology and economics, through to ecology and evolution. Their consequences are often catastrophic – for example, the sudden jamming of traffic in an information system or on a highway; a fatal change within the human immune system; so-called punctuated equilibria in evolution leading to the sudden extinction of entire species; and crashes in a financial market . The ability to generate large internal, endogenous changes is a defining characteristic of complex systems, and arguably of Nature and Life itself since it leads to evolution through innovation. Thinking through to the management of real-world complex systems and collectives, and in particular the risk associated with such catastrophic changes, one wonders whether these large events could eventually be controlled, or even avoided altogether?

Here we take the first few steps in the direction of ‘risk-management’ in collectives. We consider the specific problem of a complex adaptive multi-agent population, competing for some limited resource. The system is a variant of Arthur’s famous El Farol bar model , and exhibits large, self-generated changes. We show that information about the large endogenous changes becomes encoded in the system ahead of the large change itself. The implication is that with a reasonable amount of information, the large change would cease to be a ‘surprise’. The present work therefore has relevance to both the forward and reverse problem in collectives. In ‘pre-design’ collectives – i.e. collectives in which the individual components can be tailor-made to have certain fixed properties – such knowledge could help in the design of ‘safe’ agents whose collective behavior is such that large changes are avoided altogether. In ‘post-design’ collectives – i.e. collectives for which a designer has no control over the properties of the individual components – one could nevertheless hope to control such large changes in order to minimize the potential damage. In particular, one could imagine some form of soft monitoring, whereby an external regulator monitored the system output, detected certain precursors suggesting a dangerous build-up to a large change, then intervened to divert it. Going further, we have recently shown  that such systems can actually be ‘vaccinated’ ahead of such large changes thereby affording the system a degree of temporary immunity.

The science of complex systems, as befits its name, lacks a simple definition . It has been used to refer to systems that are intermediate between perfect order and perfect disorder; or even as a simple restatement of the cliché that the behaviour of some systems as a whole can be more than the sum of their parts. Formal definitions of complexity do exist in the computational and information sciences, but apply to specific systems . Here we focus on systems where both the properties of the individual components and the nature of their interactions are reasonably well understood. The constituents themselves can be rather simple, and the relation between any two may also be well understood, yet the collective behavior of the ensemble still manages to escape any simple explanation. A simple example which captures the idea of (co-)evolution and bounded rationality is Arthur’s El Farol bar problem . At its most basic level, this is an example of a game involving a fixed population of players (agents) $`N_{tot}`$ competing for a limited resource $`L`$. No more than $`L<N_{tot}`$ agents can win at each timestep. The agents are adaptive and try to predict the next winning outcome, which is determined only by their own choices. Each agent is equipped with a number of strategies from which to choose, and can do so adaptively based on its present environment. Recent experiments confirm that humans indeed possess such different strategies, despite being faced with the same history . In short, the El Farol bar problem contains the following key ‘complex’ elements:

* A population of heterogeneous agents with bounded rationality.
* Agents who use history-based strategies upon which they base their actions.
* A method of aggregating agent’s behaviour into a global outcome, which introduces feedback into the system.

Our study focuses on a generic multi-agent system, based on the El Farol bar problem, which has already been shown to reproduce statistical and dynamical features similar to those of a real-world complex adaptive system: a financial market. (See www.unifr.ch/econophysics for a detailed account of the extent to which financial markets represent complex systems, and their associated statistical and dynamical properties). Our model also exhibits the crucial feature of seemingly spontaneous large changes of variable duration . Although our study is therefore specific to a particular model, we believe that any multi-agent model which shares the above elements would benefit from the techniques described in this chapter. Recently the concept of an *agent* has become increasingly important in both artificial intelligence and mainstream computer science . In contrast to the agents considered within El Farol-like complex systems as in this Chapter, the agents in these other fields are typically more sophisticated and are often utilized to perform specific tasks in a wide range of different areas, e.g. e-commerce, classification, information retrieval and management of networks . By studying complex systems involving more basic agents, but which retain the ability to react to their present environment, we aim to develop approaches that may be useful in more complicated collectives comprising richer components.

The Chapter is arranged as follows. In §2 we introduce our generic complex system and in §3 demonstrate its ability to generate time-series which include occasional large movements. In §4 we investigate the game dynamics during these large changes and determine what is occurring microscopically within the model system. In §4.1 we introduce a method of understanding when a large change is possible within the system, and obtain approximate expressions for its duration in §4.2. This leads to an understanding of how such large changes may be avoided.

## 2 The Grand Canonical Bar Game

We consider a generic complex system in which a population of $`N_{tot}`$ heterogeneous agents with limited capabilities and information repeatedly compete for a limited global resource. The model is a generalization of Arthur’s bar model where agents will only take part, or play, if they are sufficiently confident of the strategies they hold. A strategy is a forecasting rule that generates a prediction of the next winning outcome based on knowledge of the recent history of winning choices. This is the Grand Canonical Bar Game (GCBG),111Grand Canonical is a term used in statistical physics to describe a system with a variable number of particles. It captures some of the key behavioural phenomena that are important in collectives/complex systems; those of competition, frustration and adaptability. It is also a ‘minimal’ system of only few parameters.

The GCBG comprises a number of agents $`N_{tot}`$ who repeatedly decide whether to enter a game where the choices are option A or B. Because of the limited global resource, only agents who are sufficiently confident of winning will participate at each timestep. The outcome at each timestep represents the winning decision, A or B. The agents are adaptive in their strategy choices, but not evolutionary; there is no discovery of new strategies by the agents. A maximum of $`L(t)`$ agents can win at each timestep. Changing $`L(t)`$ affects the system’s quasi-equilibrium; hence $`L(t)`$ can be used to mimic the changing external environment. In the limit that $`L(t)`$ is time-independent, and all agents are forced to play at each timestep, our model reverts to Arthur’s El Farol bar model. A schematic of the game structure is shown in Figure 1, clearly indicating the feedback present within the system. The method of encoding the outcomes in the game via a binary alphabet, and the associated strategy space, were introduced by Challet and Zhang . However the GCBG differs from the basic Minority Game in its use of an external resource level, and a variable number of participating agents per timestep as a result of the finite confidence level. The agents are of limited yet similar capabilities. Each agent is assigned a ‘brain-size’ $`m`$; this is the length of the past history bit-string that an agent can use when making its next decision. A common bit-string of the $`m`$ most recent outcomes is made available to the agents at each timestep. This is the only information they can use to decide which option to choose in subsequent timesteps. Each agent randomly picks $`q`$ strategies at the beginning of the game, with repetitions allowed. A strategy uses information from the historical record of winning options to generate a prediction . After each turn, the agent assigns one (virtual) point to each of his strategies which would have predicted the correct outcome, and minus one for an incorrect prediction. The resulting time-series appears ‘random’ yet is non-Markovian, with subtle temporal correlations which put it beyond any random-walk based description. Multi-agent games such as the present GCBG may be simulated on a computer, but can also be expressed in analytic form. In the next section we introduce some notation used in the remainder of this chapter.

### 2.1 Notation

A subset of the agent population, who are sufficiently confident of winning, are active at each timestep. At each timestep $`t`$ we denote the number of agents choosing option A as $`N_A(t)`$ and the number choosing B as $`N_B(t)`$. If $`L(t)N_A(t)>0`$ the winning decision is A and vice-versa. The winning decision is thus given by

$`w(t)`$ $`=`$ $`\text{sgn}\left[L(t)N_A(t)\right],`$
where $`\text{sgn}[x]`$ is the sign function defined by

$$\text{sgn}[x]=\{\begin{array}{cc}1\hfill & \text{for }x<0,\hfill \\ 0\hfill & \text{for }x=0,\hfill \\ 1\hfill & \text{for }x>0.\hfill \end{array}$$

We denote by $`w(t)`$ the winning option at time $`t`$, where a value of $`1`$ option A, and $`1`$ option B. We frequently have to represent the two possible choices, A or B, in numerical form, and we use the encoding that $`1`$ always implies option A and $`1`$ implies option B. If $`w(t)=0`$, indicating no clear winning option, this value is replaced with a random coin toss.

The ‘excess demand’ $`D(t)=N_A(t)N_B(t)`$ (which mimics price-change in a market) and number $`V(t)=N_A(t)+N_B(t)`$ of active agents (which mimics volume) represent output variables. These two quantities fluctuate with time, and can be combined to construct other global quantities of interest for the complex system studied. We define the cumulative excess demand as $`P(t)`$. In the context of a financial market, this can be regarded as a pseudo-price. Typically we use the example of a financial market, but the excess demand can be interpreteted in many different circumstances, e.g. as a measure of resource utilization within a system. If we define $`L(t)=\varphi V(t)`$, where $`0\varphi 1`$, then only a fraction $`\varphi `$ of the active population can win. By varying the value of $`\varphi `$, it is possible to change the equilibrium value of $`D(t)`$, see Figure 2.

The only global information available to the agents is a common bit-string ‘memory’ of the $`m`$ most recent outcomes. Consider $`m=2`$; the $`P=2^m=4`$ possible history bit-strings are $`AA`$, $`AB`$, $`BA`$ and $`BB`$, with the right-most letter representing the winning choice at the last timestep. The history, or equivalently, the global information available to the agents, can be represented in decimal form $`\mu \{0,1,\mathrm{},P1\}`$:

$`\mu (t)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{m}{}}}2^{i1}\left[w(ti)+1\right].`$
The global information $`\mu (t)`$ is updated by dropping the first bit and concatenating the latest outcome to the end of the string.

At the beginning of the game, each agent randomly picks $`q`$ ($`>1`$) strategies, making the agents heterogeneous in their strategy sets.222Agents are only adaptive if they have more than one strategy to play with; for $`q=1`$ the game has a trivial periodic structure. This initial strategy assignment is fixed from the outset of each simulation and provides a systematic disorder which is built into each run. This is referred to as the *quenched disorder* present in the game. An agent decides which option to choose at a given time step based on the prediction of a strategy, which consists of a response, $`a^{\mu (t)}\{1,1\}`$ to the global information $`\mu (t)`$. For its current play, an agent chooses its strategy that would have performed best over the history of the game until that time, i.e. has the most virtual points.

Agents have a time horizon $`T`$ over which virtual points are collected and a threshold probability level $`\tau `$ which mimics a ‘confidence’. Only strategies having $`r`$ points are used, where

$`r`$ $`=`$ $`T(2\tau 1).`$
We call these *active* strategies. Agents with no active strategies within their individual set of $`q`$ strategies do not play at that timestep and become temporarily inactive. Agents with one or more active strategies play the one with the highest virtual point score; any ties between active strategies are resolved using a coin-toss. If an agent’s threshold to play is low, we would expect the agent to play a large proportion of the time as their best strategy will have invariably scored higher than this threshold. Conversely, for high $`\tau `$, the agent will hardly play at all. The coin-tosses used to resolve ties in decisions (i.e. $`N_A=N_B`$) and active-strategy scores, inject stochasticity into the game’s evolution. The implementation of the strategies is discussed in §2.2.

After each turn, agents update the scores of their strategies with the reward function

$`\chi [N_A(t),L(t)]=\text{sgn}\left[N_A(t)L(t)\right],`$ (1)
namely $`+1`$ for choosing the correct/winning outcome, and $`1`$ for choosing the incorrect outcome. The virtual points for strategy $`R`$ are updated via

$`S_R(t)`$ $`=`$ $`{\displaystyle \underset{i=tT}{\overset{t1}{}}}a_R^{\mu (i)}\chi [N_A(i),L(i)],`$
where $`a_R^{\mu (t)}`$ is the response of strategy $`R`$ to global information $`\mu (t)`$, and the summation is taken over a rolling window of fixed length $`T`$.333It is also possible to construct an exponentially weighted window of characteristic length $`T`$, with a decay parameter $`\lambda =11/T`$. The strategy score updating equation (4) becomes    $`𝐒(t+1)`$ $`=`$ $`𝐚^{\mu (t)}w(t)+\lambda 𝐒(t).`$   This has the advantage of removing the hard cut off from the rolling window (a Fourier transform of the demand can show the effect of this periodicity in the time series), and is a rapid recursive calculation. But strategy scores are no longer integers, and the probability of a strategy score tie is now very low. Since one source of stochasticity has now been effectively removed, there must be occasional periods of inactivity to inject stochasticity which helps to prevent group, or cyclic, behaviour when the game traces out a deterministic path. To start the simulation, we set the initial strategy scores to be zero: $`S_R(0)=0`$. Because of the feedback in the game, any particular strategy’s success is short-lived. If all the agents begin to use similar strategies, and hence make the same decision, such a strategy ceases to be profitable and is subsequently dropped. This encourages heterogeneity amongst the agent population.

The demand $`D(t)`$ and volume $`V(t)`$, which can be identified as the output from the model system, are given by

$`D(t)`$ $`=`$ $`𝐧(t)\left(𝐚^{\mu (𝐭)}\left[𝐒(t)r\right]\right)={\displaystyle \underset{R=1}{\overset{Q}{}}}n_R(t)a_R^{\mu (t)}\left[S_R(t)r\right],`$ (2a)
$`V(t)`$ $`=`$ $`𝐧(t)\left[𝐒(t)r\right]={\displaystyle \underset{R=1}{\overset{Q}{}}}n_R(t)\left[S_R(t)r\right].`$ (2b)
where $`n_R(t)`$ represents the number of agents playing strategy $`R`$ at timestep $`t`$.

The demand $`D(t)`$ is made up of two groups of agents at each timestep: $`D_D(t)`$ agents who act in a deterministic manner, i.e. do not require a coin toss to decide which choice to make - this is because they either have one strategy that is better than their others, or because their highest-scoring strategies are tied but give the same response to the history $`\mu (t)`$, and $`D_U(t)`$ agents that act in an ‘undecided’ way, i.e. they require the toss of a coin to decide which choice to make - this is because they have two (or more) highest-scoring tied strategies which give different responses at that turn of the game. Inactive agents do not contribute to the demand. Hence we can rewrite (2a) as

$`D(t)`$ $`=`$ $`D_D(t)+D_U(t).`$ (3)
Without the stochastic influence of the undecided agents the game will tend to exhibit group, or cyclic behaviour, where the game eventually traces out a deterministic path. The period of this cyclic behaviour is dependent on the quenched disorder present in the simulation, and could be very long.

The number of agents holding a particular combination of strategies can also be expressed as a $`q`$-dimensional tensor $`𝛀`$ , where the entry $`\mathrm{\Omega }_{R1,R2,\mathrm{}}`$ represents the number of agents holding strategies $`\{R1,R2,\mathrm{}\}`$. This quenched disorder is fixed at the beginning of the game. It is useful to construct a symmetric configuration $`𝚿`$ in the sense that $`\mathrm{\Psi }_{R1,R2,\mathrm{}}=\mathrm{\Psi }_{p\{R1,R2,\mathrm{}\}}`$ where $`p\{R1,R2,\mathrm{}\}`$ is any permutation of the strategies $`R1,R2,\mathrm{}`$; for $`q=2`$ we let $`𝚿=\frac{1}{2}(𝛀+𝛀^T)`$ . Elements $`\mathrm{\Psi }_{R,R^{}}`$ enumerate the number of agents holding both strategy $`R`$ and $`R^{}`$. We focus on $`q=2`$ strategies per agent, although the formalism can be generalised. At timestep $`t`$, $`D_D(t)`$ can now be expressed as

$`D_D(t)`$ $`=`$ $`{\displaystyle \underset{R=1}{\overset{Q}{}}}a_R^{\mu (t)}[S_R(t)r]{\displaystyle \underset{R^{}=1}{\overset{Q}{}}}\left(1+\text{sgn}\left[S_R(t)S_R^{}(t)\right]\right)\mathrm{\Psi }_{R,R^{}}.`$
The number of undecided agents $`N_U`$ is given by

$`N_U(t)`$ $`=`$ $`{\displaystyle \underset{R,R^{}}{}}[S_R(t)r]\delta (S_R(t)S_R^{}(t))[1\delta (a_R^{\mu (t)}a_R^{}^{\mu (t)})]\mathrm{\Psi }_{R,R^{}}`$
and hence the demand of the undecided agents $`D_U(t)`$ is distributed binomially:

$`D_U(t)`$ $`=`$ $`2\text{Bin}(N_U(t),{\displaystyle \frac{1}{2}})N_U(t)`$
where $`\text{Bin}(n,p)`$ is a sample from a binomial distribution of $`n`$ trials with probability of success $`p`$.

### 2.2 The strategy space

The strategy space analysis in the Section was inspired by the work of Challet and Zhang . A strategy consists of a response, $`a^\mu \{1,1\}`$ to each possible bit-string $`\mu `$, $`a^\mu =1`$ option A, and $`a^\mu =1`$ option B. Consider $`m=2`$, each strategy can be represented by a string of $`P=4`$ bits $`[\text{i j k l}]`$ with $`i,j,k,l=1`$ or $`+1`$ corresponding to the decisions based on the histories AA, AB, BA and BB respectively. For example, strategy $`[\text{-1 -1 -1 -1}]`$ corresponds to deciding to pick option B irrespective of the $`m=2`$ bit-string. $`[\text{1 1 1 1}]`$ corresponds to deciding to pick option A irrespective of the $`m=2`$ bitstring. $`[\text{1 -1 1 -1}]`$ corresponds to deciding to pick option A given the histories AA or BA, and pick option B given the histories AB or BB. A subset of strategies can further be classed as one of the following:

* *anti-correlated*: for example, any two agents using the strategies $`[\text{-1 -1 -1 -1}]`$ and $`[\text{1 1 1 1}]`$ would take the opposite action irrespective of the sequence of previous outcomes. Hence one agent will always do the opposite of the other agent, and their net effect on the excess demand $`D(t)`$ will be zero.
* *un-correlated*: for example, any two agents using the strategies $`[\text{-1 -1 -1 -1}]`$ and $`[\text{1 -1 1 -1}]`$ would take the opposite action for two of the four histories, and take the same action for the remaining two histories. Assuming that the $`m=2`$ histories occur equally often, the actions of the two agents will be uncorrelated on average

A convenient measure of the distance of any two strategies is the relative Hamming distance, defined as the number of bits that need to be changed in going from one strategy to another. For example, the Hamming distance between $`[\text{-1 -1 -1 -1}]`$ and $`[\text{1 1 1 1}]`$ is 4, while the Hamming distance between $`[\text{-1 -1 -1 -1}]`$ and $`[\text{1 -1 1 -1}]`$ is just 2.

The collection of all possible strategies (and their associated virtual points) is hereafter referred to as the *strategy space*, see e.g. Figure 3. This object can be thought of as a common property of the game itself, being updated centrally after each timestep. Each individual agent monitors only a fixed subset $`q`$ of the $`Q`$ possible strategies (with the small caveat that a repeated strategy choice is possible). After each turn, an agent assigns one (virtual) point to each of its strategies which would have predicted the correct outcome. The virtual points for each strategy can be represented by the strategy score vector $`𝐒(t)`$, and are given by

$`𝐒(t)`$ $`=`$ $`{\displaystyle \underset{i=tT}{\overset{t1}{}}}𝐚^{\mu (i)}\chi [N_A(i),L(i)]`$ (4)
In total, there are $`Q=2^P`$ possible strategies which define the decisions in response to all possible $`m`$ history bit-strings. This is referred to as the full strategy space (FSS). However, the principal features of the system are reproduced in a smaller reduced strategy space (RSS) of $`Q=2^{m+1}`$ strategies wherein any two strategies are either un-correlated or anti-correlated , i.e. separated by a Hamming distance of either $`2^m`$ or $`2^{m1}`$. The full and reduced strategy space $`𝐚`$ for $`m=2`$ has been reproduced in Figure 4. Each row represents a strategy, each column is assigned to a particular history, giving the strategy space a dimension of $`P\times Q`$. The prediction of strategy $`R`$ to information $`\mu `$ is $`a_R^\mu `$ and corresponds to the $`(R,\mu )`$ element of $`𝐚`$. The ordering of the rows in unimportant.

Within the strategy space, each strategy $`R`$ has an anticorrelated strategy $`\overline{R}`$. We note that the anticorrelated strategies are effectively redundant, as their predictions and strategy scores can be recovered from their anticorrelated pair:

$`𝐚_{\overline{R}}^\mu `$ $`=`$ $`𝐚_R^\mu ,`$ (5a)
$`𝐒_{\overline{R}}(t)`$ $`=`$ $`𝐒_R(t).`$ (5b)
Equation (5a) is true by definition, i.e. a strategy and its anticorrelated pair always give the opposite prediction. However (5b) requires a symmetric scoring rule to be used, which is satisfied by (1). Thus we can reproduce the dynamics using a space of just $`P`$ strategies, in which agents choose both a strategy and whether to agree or disagree with its prediction. Reducing the size of the strategy space is advantageous as it reduces memory requirements and increases the speed of the simulation. For a more detailed description of how to implement such a system, see .

## 3 Demonstration of large changes

A ubiquitous feature of complex systems is that large changes, or ‘extreme events’, arise far more often than would be expected if the individual agents acted independently. We frequently refer to crashes, but are interested in large moves in either direction. With a suitable choice of parameters the GCBG is able to generate time series which include occasional large movements, see e.g. Figure 5. The game can be broadly classified into three regimes:

1. The number of strategies in play is much greater than the total available: groups of traders will play using the same strategy and therefore crowds should dominate the game .
2. The number of strategies in play is much less than the total available: grouping behaviour is therefore minimal.
3. The number of strategies in play is comparable to the total number available.

We focus on the third regime, since this yields seemingly random dynamics with occasional large movements.

Large changes seem to exhibit a wide range of possible durations and magnitudes making them difficult to capture using traditional statistical techniques based on one or two-point probability distributions. A common feature, however, is an obvious trend (i.e. to the eye) in one direction over a reasonably short time window: we use this as a working definition of a large change. (In fact, all the large changes discussed here represent $`>3\sigma `$ events.) In both our model and the real-world system, these large changes arise more frequently than would be expected from a random-walk model .

In §3.1 we consider the distribution of the excess demand created by our generic system, and perform a simple analysis to discuss the statistics of extreme events. To determine whether large movements occur due to a single random event, we examine the stochastic influence present within the model in §3.2 and find that they arise through a global cooperation occurring over the whole system.

### 3.1 Tail estimation

Traditional parametric statistical methods are ill-suited for dealing with extreme events, which have little historical data associated with them. Provided the distribution has a finite variance and the increments are independent, the Central Limit Theorem will apply near the centre of the distribution, but does not tell us anything about the tails .444The speed at which the distribution will converge to a Gaussian is given by the Berry-Esseen theorem; a finite second moment ensures Gaussian behaviour, and the speed of convergence is controlled by the size of the third moment of $`|x|`$ . In recent years extreme value theory (EVT), a branch of probability theory which focuses explicitly on extreme outcomes, has received increasing attention . EVT considers the distribution of extreme returns rather than the distribution of all returns. It can be shown that the limiting distribution of extreme returns observed over a long time-period is largely independent of the distribution of returns itself. The upper tail of a fat-tailed cumulative distribution function $`F`$ behaves asymptotically like the tail of the Pareto distribution given by $`1F(x)cx^\alpha `$, for $`c>0`$, $`\alpha >0`$ and $`xC`$, where $`C`$ is a threshold above which the assumed algebraic form is valid, and $`c`$ is a normalising constant. The tail index $`\alpha `$ determines the heaviness of the tail of a distribution and plays a key role in tail-related risk measures, representing the maximal order of finite moments. Only the first $`k`$ moments, where $`k<\alpha `$, are bounded. The greater the tail index, the ‘fatter’ the tail, and the greater the incidence of extreme events.

An important issue in the study of fat-tailed distributions is the estimation of the tail index $`\alpha `$. There are a number of methods to estimate the tail index, some using asymptotic results from EVT, from which values can be estimated using maximum likelihood techniques. However the Hill estimator is commonly used, as it is suitable for tail estimation of fat-tailed distributions and is relatively easy to implement . Consider a sample $`X_1,X_2,\mathrm{},X_n`$ of $`n`$ observations drawn from a stationary iid process, and let $`X_{(1)}X_{(2)}\mathrm{}X_{(n)}`$ be the descending order statistics. The Hill estimator is based on the difference between the $`m`$th largest observation and the average of the $`m`$ largest observations:

$`\widehat{\xi }_m`$ $`=`$ $`{\displaystyle \frac{1}{m}}{\displaystyle \underset{i=1}{\overset{m}{}}}\mathrm{log}X_{(i)}\mathrm{log}X_{(m)}`$
where $`\xi =1/\alpha `$ is the shape parameter and $`m`$ is the number of order statistics used in the tail estimation. The appropriate choice of value $`m`$ is a non-trivial problem, since this requires us to decide where the tail begins. There is a tradeoff between the bias and variance of the estimator in choosing $`m`$. If we choose a large value of $`m`$, the number of order statistics used increases and the variance of the estimator will decrease. However, choosing a high $`m`$ also introduces some observations from the centre of the distribution and the estimation becomes biased. But if it is too small, the estimate will be based on a just a few of the largest observations and the estimator will lack precision. Several methods for the determination of an optimal sample fraction for the Hill estimator have been proposed . In Figure 6 a Hill plot is constructed from the simulation data used to create Fig. 5. A threshold is selected from the plot where the shape parameter is fairly stable, giving an estimate of the tail index $`\alpha 3.7`$. To obtain the most accurate tail estimate, a combination of several techniques should be considered.

### 3.2 The effect of stochasticity

To investigate the stochastic influence due to the effect of coin tosses, we have studied their occurrence during the simulation. We define two types of coin toss: type (i), which occurs when an agent has a tie between active strategy options which predict differing outcomes, and type (ii), which occurs when the number of agents choosing option A is equal to the number choosing option B. The main conclusion is that there is no single-coin toss that immediately causes a large change within the system. and The large movements observed as arising endogenously in the system, result from the organisation of temporal and spatial (i.e. strategy space) correlations. In common with other complex adaptive systems, this organisation does not arise in general from a nucleation phase diffusing across the system, e.g. it cannot be traced to a particular coin-toss by a particular agent which triggers the “avalanche”. Rather, it results from a progressive and more global cooperation occurring over the whole system via repetitive interactions .

## 4 A study of extreme events

In this section we consider the generic complex system introduced in §2 in which a population of $`N_{tot}`$ heterogeneous agents with limited capabilities and information, repeatedly compete for a limited global resource. We focus on extreme events which are endogenous (i.e. internally-produced by the agents themselves) and provide a microscopic understanding as to their build-up and likely duration.

### 4.1 Nodal weight decomposition

The dynamics of the game history can be usefully represented on a directed de Bruijn graph . This is a graph whose nodes are sequences of symbols from some alphabet, and whose edges represent possible transitions between these nodes. It has numerous interesting properties and is frequently discussed in the context of parallel algorithms and communications networks. This is an effective method of representing the evolution of the system, and in Figure 7 we plot the de Bruijn graph for an $`m=3`$ game.

Large changes occur when connected nodes become persistent and the game makes successive moves in the same direction. Only nodes 0 and 7 can exhibit perfect nodal persistence, where an allowed transition can return the system to exactly the same node. This is the simplest type of large change, e.g. $`\mu (t)=0,0,0,0,\mathrm{}`$, where all successive price changes are in the *same* direction. We call this a ‘fixed-node crash’ (or rally). There are many other possibilities reflecting the wide range of forms and durations that a large change can undertake. For example, on the $`m=3`$ de Bruijn graph in Fig. 7 the cycle $`\mu =0,0,1,2,4,0,\mathrm{}`$ has four out of the five transitions producing price-changes of the same sign (it is persistent on nodes 1, 2, 4 and antipersistent on node 0). We call this a ‘cyclic-node crash’. Stable behaviour occurs on a path where all transitions are equally visited, e.g. $`\mu (t)=0,0,1,3,6,5,3,7,7,6,4,1,2,5,4,\mathrm{}`$.

To identify moments when large changes occur, we need to recognise when history nodes are likely to become persistent. Whether a node is persistant or not will depend on the action of the agents at that timestep, which is determined by the predictions of the strategies they hold. If the majority of active strategies generate a prediction of the same outcome, then the demand at that timestep is likely to continue in that direction. Thus a suitable condition for a large movement is when there is a concensus of opinion regarding the next outcome amongst the active strategies. This occurs when the pattern of active strategies within the strategy space matches up with the pattern of stratgy predictions for a history node $`\mu `$, see Figure 8. The strategy predictions will depend on the global information $`\mu (t)`$, and the distribution of active strategies depends on the strategy score vector $`𝐒(t)`$. At each timestep $`𝐒(t)`$ is updated according to whether a strategy predicted the winning outcome $`w(t)`$; the incremental strategy score is given by a column of the strategy matrix $`𝐚`$ determined by $`\mu (t)`$. In total there are $`P`$ orthogonal increment vectors $`𝐚^\mu `$, one for each value of $`\mu `$. We can express the strategy score (4) as

$`𝐒(t)=c_0𝐚^0+c_1𝐚^1+\mathrm{}+c_{P1}𝐚^{P1}={\displaystyle \underset{j=0}{\overset{P1}{}}}c_j𝐚^j,`$ (6)
where $`c_j`$ represents the *nodal weight* for history node $`\mu =j`$. The nodal weights represent the number of negative return transitions from node $`\mu `$ minus the number of positive return transitions, in the time window $`tTt1`$. The values of these nodal weights are important in identifying periods in which large changes can occur. If the value of a nodal weight is near zero, this implies that the number of active strategies predicting outcome A will be similar to the number predicting outcome B when that node is reached. The excess demand will then depend on the quenched disorder and the undecided agents, but is likely to be small. Conversely, if the value of a nodal weight is significantly different from zero, this indicates a bias in the strategy space, i.e. the same pattern is evident in the active strategies and the strategy predictions. When the game trajectory hits this node, the majority of active strategies will predict the same outcome, and there will be a consensus amongst the agents. This is likely to lead to a large excess demand, and therefore a large change. Thus high absolute nodal weight implies persistence in transitions from that node, i.e. persistence in $`D|\mu `$.

The mean of the strategy scores predicting an $`A`$ or $`B`$ at node $`\mu `$ are linked to their nodal weight:

$`\text{E}\left[S_{Ra_R^\mu =1}(t)\right]`$ $`=`$ $`c_\mu (t),`$ (7a)
$`\text{E}\left[S_{Ra_R^\mu =1}(t)\right]`$ $`=`$ $`c_\mu (t).`$ (7b)
This symmetry occurs because each strategy has an anticorrelated strategy present in the strategy space.

The nodal weight decomposition provides a succinct method of describing the state of the strategy space. We are concerned with high nodal weights on game cycles that can become persistent. When a nodal weight value becomes large, this is a warning sign that the system is in a suitable state to undergo a large change. Figure 9 illustrates a large change which starts as a fixed-node crash, and then subsequently becomes a cyclic-node crash.

### 4.2 Estimating the crash length

We are interested in the dynamics of large changes, and use a simplified version of the system to obtain an analytic expression for the expected crash length. Reference  showed that the GCBG can be usefully described as a stochastically disturbed deterministic system. As stated in (3), the demand can be divided into contributions arising from agents acting on a definite strategy $`D_D(t)`$, and that from the undecided agents $`D_U`$. The average contribution of the undecided agents to the net demand will be zero, i.e. $`\text{E}[D_U(t)]=0`$. Averaging over our model’s stochasticity in this way yields a description of the game’s deterministic dynamics. By examining these dynamics, we can determine when a large-change is likely. For $`q=2`$ the demand function can be expressed as

$`D(t)`$ $`=`$ $`{\displaystyle \underset{R=1}{\overset{Q}{}}}a_R^{\mu (t)}[S_R(t)r]{\displaystyle \underset{R^{}=1}{\overset{Q}{}}}\left(1+\text{sgn}\left[S_R(t)S_R^{}(t)\right]\right)\mathrm{\Psi }_{R,R^{}},`$ (8)
where $`𝚿`$ is the symmetrized strategy allocation matrix which constitutes the quenched disorder present during the system’s evolution. The volume $`V(t)`$ is given by the same expression as $`D(t)`$ replacing $`a_R^{\mu (t)}`$ by unity.

For the parameter ranges of interest, the choice about whether a strategy is played by an agent is more determined by whether that strategy’s score is above the threshold, than whether it is their highest-scoring strategy. This is because agents are only likely to have at most one strategy whose score lies above the threshold for confidence levels $`r0`$. Making the additional numerically-justified approximation of small quenched disorder (i.e. the variance of the entries in the strategy allocation matrix $`𝚿`$ is smaller than their mean for the parameter range of interest ), the demand and volume become

$`D(t)`$ $`=`$ $`{\displaystyle \underset{R=1}{\overset{Q}{}}}a_R^{\mu (t)}[S_R(t)r]{\displaystyle \underset{R^{}=1}{\overset{Q}{}}}{\displaystyle \frac{N}{Q^2}}`$ (9a)
$`=`$ $`{\displaystyle \frac{N}{Q}}{\displaystyle \underset{R=1}{\overset{Q}{}}}a_R^{\mu (t)}\left[S_R(t)r\right]`$
$`=`$ $`{\displaystyle \frac{N}{Q}}{\displaystyle \underset{R=1}{\overset{Q}{}}}a_R^{\mu (t)}{\displaystyle \frac{1}{2}}\left(1+\text{sgn}\left[S_R(t)r\right]\right)`$
$`=`$ $`{\displaystyle \frac{N}{2Q}}{\displaystyle \underset{R=1}{\overset{Q}{}}}a_R^{\mu (t)}\text{sgn}\left[S_R(t)r\right],`$
$`V(t)`$ $`=`$ $`{\displaystyle \frac{N}{Q}}{\displaystyle \underset{R=1}{\overset{Q}{}}}\left[S_R(t)r\right]={\displaystyle \frac{N}{2}}+{\displaystyle \frac{N}{2Q}}{\displaystyle \underset{R=1}{\overset{Q}{}}}\text{sgn}\left[S_R(t)r\right].`$ (9b)
Let us suppose persistence on node $`\mu =0`$ starts at time $`t_0`$. How long will the resulting large change last? To answer this, we decompose (9a) into strategies which predict option $`A`$ at $`\mu =0`$, and those that predict $`B`$. We first consider the particular case where the node $`\mu =0`$ was not visited during the previous $`T`$ timesteps, hence the loss of score increment from time-step $`tT`$ will not affect $`𝐒(t)`$ on average. At any later time $`t_0+\tau `$ during the large change, (9a) and (9b) are given by

$`D(t_0+\tau )`$ $`=`$ $`{\displaystyle \frac{N}{2Q}}\left\{{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R(t_0)r\tau \right]{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R(t_0)r+\tau \right]\right\},`$ (10)
$`V(t_0+\tau )`$ $`=`$ $`{\displaystyle \frac{N}{2}}+{\displaystyle \frac{N}{2Q}}\left\{{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R(t_0)r\tau \right]+{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R(t_0)r+\tau \right]\right\}.`$
The magnitude of the demand $`|D(t_0+\tau )|`$ decreases as the persistence time $`\tau `$ increases, and the large change will end at time $`t_0+\tau _c`$ when the right-hand side of (10) becomes zero. We denote the persistence time or ‘crash-length’ by $`\tau _c`$. It is easy to obtain an upper limit for the duration of a large change by determining when the demand $`D`$ changes sign. This occurs when the mean of the scores of the strategies predicting $`A`$ returns to 0, i.e.

$`\tau _c`$ $`=`$ $`\text{E}\left[S_{Ra_R^\mu =1}(t_0)\right]`$
$`=`$ $`c_0(t_0)`$
from (7b). Thus the crash length will depend on the score difference between strategies predicting option $`A`$ and those predicting $`B`$. In the more general case, where the node $`\mu =0`$ was visited during the previous $`T`$ timesteps, $`\tau _c`$ is given by the largest $`\tau `$ value which satisfies

$`\tau `$ $`=`$ $`\left(c_0(t_0)+{\displaystyle \underset{\{t^{}\}}{}}sgn\left[D(t^{})\right]\right)`$ (11)
where $`\{t^{}\}(\mu (t^{})=0t_0Tt^{}t_0T+\tau )`$. The summation accounts for any $`\mu =0`$ transitions in the period $`t_0T`$ to $`t_0T+\tau `$.

To obtain an improved estimate of $`\tau _c`$, and investigate the behaviour of the demand and volume, it is necessary to assume a distribution for the strategy scores. In Figure 10 we plot an example of the strategy score distribution prior to a large change. We assume that the scores have a near-Normal distribution, i.e. $`S_{Ra_R^\mu =1}(t_0)N[c_0(t_0),\sigma ]`$ as shown Fig. 11.555At some points during the simulation the distribution can be more complex, e.g. multimodal behaviour. Consequently, prior to a large change, the score distribution tends to split into two halves. Substitution into (10) gives the expected demand and volume during the crash:

$`D(t_0+\tau )`$ $``$ $`\mathrm{\Phi }\left[{\displaystyle \frac{c_0(t_0)r+\tau }{\sigma }}\right]\mathrm{\Phi }\left[{\displaystyle \frac{c_0(t_0)r\tau }{\sigma }}\right],`$
$`V(t_0+\tau )`$ $``$ $`\mathrm{\Phi }\left[{\displaystyle \frac{c_0(t_0)r+\tau }{\sigma }}\right]+\mathrm{\Phi }\left[{\displaystyle \frac{c_0(t_0)r\tau }{\sigma }}\right].`$
These forms are illustrated in Figure 12. As the spread in the strategy score distribution is increased, the dependence of $`D`$ and $`V`$ on the parameters $`\tau `$ and $`r`$ becomes weaker and the surfaces flatten out leading to a smoother drawdown, as opposed to a sudden severe crash. In the limit $`\sigma 0`$, the crash-length $`\tau _c=c_0(t_0)|r|`$. As the parameters $`c_0,\sigma ,r`$ are varied, it can be seen that the behaviour of the demand and volume during the large change can exhibit markedly different qualitative forms.

### 4.3 Repeated occurence of large changes

We now turn to the important practical question of whether history will repeat itself, i.e. given that a large change has recently happened, is it likely to happen again? If so, is it likely to be even bigger? Suppose the system has built up a negative nodal weight for $`\mu =0`$ at some point in the game (see Fig. 13a). It then hits node $`\mu =0`$ at time $`t_0`$ producing a large change (Fig. 13b). The nodal weight $`c_0`$ is hence restored to zero (Fig. 13c). In this model the previous build-up is then forgotten because of the finite $`T`$ score window, hence $`c_0`$ becomes positive (Fig. 13d). The system then corrects this imbalance (Fig. 13e), restoring $`c_0`$ to 0. The large change is then forgotten, hence $`c_0`$ becomes negative (Fig. 13f). The system should therefore crash again - however, a crash will *only* re-appear if the system’s trajectory subsequently returns to node $`\mu =0`$. Interestingly, we find that the *disorder* in the initial distribution of strategies among agents (i.e. the quenched disorder in $`𝚿`$) can play a deciding role in the issue of ‘births and revivals’ of large changes since it leads to a slight bias in the outcome, and hence the subsequent transition, at each node. In certain configurations, this system may demonstrate repeated instabilities leading to large changes. When $`c_{\mu (t)}=0`$ (see Fig. 13c), it follows that $`sgn[D(t)]`$ is more likely to be equal to $`sgn\left[𝐚^{\mu (t)}𝐱\right]`$ where $`x_R=_R^{}\mathrm{\Psi }_{R,R^{}}`$ is a strategy weight vector with $`x_R`$ corresponding to the number of agents who hold strategy $`R`$ . The quenched disorder therefore provides a crucial bias for determining the future trajectory on the de Bruijn graph when the nodal weight is small, and hence can decide whether a given large change recurs or simply disappears. The quenched disorder also provides a *catalyst* for building up a very large change .

## 5 Concluding remarks

We have addressed the issue of understanding, predicting and eventually controlling catastrophic endogenous changes in a collective. By utilizing information about the strategy weights within the model, we have developed a method for determining when a large change is likely within a generic complex system, the so-called GCBG, and obtained an analytic expression for its expected duration. Our work opens up the study of how a ‘complex-systems-manager’ might use this information to control the long-term evolution of a complex system by introducing, or manipulating, such large changes . As an example, we give a quick three-step solution to prevent large changes: (1) use the past history of outcomes to build up an estimate of the score vector $`\text{S}(t)`$ and the nodal weights $`\{c_{\mu (t)}\}`$ on the various critical nodes, such as $`\mu =0`$ in the case of the fixed-node crash. (2) Monitor these weights to check for any large build-up. (3) If such a build-up occurs, step in to prevent the system hitting that node until the weights have decreased. Finally we note that it can sometimes be beneficial to induce small changes ahead of time, in order to avoid larger changes in the future. We call this process ‘immunization’ and refer to Ref.  for more details.

We are grateful to David Wolpert, Kagan Tumer and Damien Challet for many useful discussions about Collectives.


text:
# Dynamical Instabilities in Extrasolar Planetary Systems

## 1. Background

The known extrasolar planets (Fig. 1 left) can be roughly divided into two groups: those with short-period, nearly circular orbits ($`a<0.07`$AU) and those with wider and more eccentric orbits ($`a>0.07`$AU). Many of the short-period planets, like their prototype 51 Peg, are so close to their parent star that tidal dissipation would have likely circularized their orbits, even if they were originally eccentric (Rasio et al. 1996). Thus, their small observed eccentricities do not provide a good indicator of their dynamical history. However, the large eccentricities of most planets with longer periods require an explanation. A planet that would have formed from a protoplanetary disk in the standard manner is unlikely to have developed such a large eccentricity, since dissipation in the disk tends to circularize orbits. Dynamical instabilities leading to the ejection of one planet while retaining another planet of comparable mass could explain the observed distribution of eccentricities.

In previous papers (Rasio & Ford 1996; Ford, Havlickova, & Rasio 2001), we conducted Monte Carlo integrations of planetary systems containing two equal mass planets initially placed just inside the Hill stability limit (Gladman 1993). We numerically integrated the orbits of such systems until there was a collision, or one planet was ejected from the system, or some maximum integration time was reached. The two most common outcomes were collisions between the two planets, producing a more massive planet in a nearly circular orbit between the two initial orbits, and ejections of one planet from the system while the other planet remains in a tighter orbit with a large eccentricity. The relative frequency of these two outcomes depends on the ratio of the planet radius to the initial semi-major axis. While this model could naturally explain how planets might acquire large eccentricities, upon comparing our results to the observed planets we found two important differences. First, for the relevant radii and semi-major axes, we found that collisions are more frequent in our simulations than nearly circular orbits among the presently known extrasolar planets. Second, of the systems which ended with one planet being ejected from the system, the eccentricity distribution of the remaining planet was concentrated in a narrow range which is greater than the typical eccentricity of the known extrasolar planets (See Fig. 1, right). Here we present results for new simulations involving two planets with unequal masses and three planets with equal masses.

## 2. Two Planets, Unequal Masses

### 2.1. Numerical Setup

The new simulations use the mixed variable symplectic algorithm (Wisdom & Holman 1991) modified to allow for close encounters between planets as implemented in the publicly available code Mercury (Chambers 1999). The results presented below are based on $`10^4`$ numerical integrations.

Throughout the integrations, close encounters between any two bodies were logged, allowing us to present results for any values of the planetary radii using a single set of orbital integrations. We consider a range of radii to allow for the uncertainty in both the physical radius and the effective collision radius allowing for dissipation in the planets. When two planets collided, mass and momentum conservation were used to compute the final orbit of the resulting single planet.

Each run was terminated when one of the following four conditions was encountered: (i) one of the two planets became unbound (which we defined as having a radial distance from the star of $`1000a_{1,\mathrm{init}}`$); (ii) a collision between the two planets occurred assuming $`R_i/a_{1,\mathrm{init}}=R_{\mathrm{min}}/a_{1,\mathrm{init}}=0.1R_{\mathrm{Jup}}/5\mathrm{AU}=0.95\times 10^5`$; (iii) a close encounter occurred between a planet and the star (defined by having a planet come within $`r_{\mathrm{min}}/a_{1,\mathrm{init}}=10R_{}/1\mathrm{AU}=0.06`$ of the star); (iv) the integration time reached $`t_{\mathrm{max}}=510^6210^7`$ depending on the masses of the planets. These four types will be referred to as “collisions,” meaning a collision between the two planets, “ejections,” meaning that one planet was ejected to infinity, “star grazers,” meaning that one planet had a close pericenter passage, and “two planets,” meaning that two bound planets remained in a (possibly new) dynamically stable configuration.

Our numerical integrations were performed for a system containing two planets, with mass ratios $`10^4<m_i/M<10^2`$, where $`m_i`$ is the mass of one of the planets and $`M`$ is the mass of the central star. A mass ratio of $`m/M10^3`$ corresponds to $`m1M_{\mathrm{Jup}}`$ for $`M=1M_{}`$. The initial semimajor axis of the outer planet ($`a_{2,\mathrm{init}}`$) was drawn from a uniform distribution ranging from $`0.9a_{1,\mathrm{init}}\left(1+\mathrm{\Delta }_c\right)`$ to $`a_{1,\mathrm{init}}\left(1+\mathrm{\Delta }_c\right)`$, where $`1+\mathrm{\Delta }_c`$ is the critical ratio above which Hill stability is guaranteed for initially circular coplanar orbits (Gladman 1993). The initial eccentricities were distributed uniformly in the range from 0 to 0.05, and the initial relative inclination in the range from $`0^{}`$ to $`2^{}`$. All remaining angles (longitudes and phases) were randomly chosen between 0 and $`2\pi `$. Throughout this paper we quote numerical results in units such that $`G=a_{1,\mathrm{init}}=M=1`$. In these units, the initial orbital period of the inner planet is $`P_12\pi `$.

### 2.2. Results

#### Branching Ratios

For the range of $`m_i/M`$ considered (up to $`10^2`$), the frequency of the various outcomes is relatively insensitive to the ratio of the total planet mass to the stellar mass. However the frequencies of the different outcomes are significantly effected by the planetary mass ratio and by the planetary radii. Collisions are most common for large planetary radii and mass ratios near unity, while ejections are most common for small planetary radii and mass ratios far from unity. The frequency of outcomes does not depend on whether the more massive planet is initially in the inner or outer orbit. In Fig. 2 (left), we show the frequency of ejection or collision as a function of the planetary radii with solid and dashed lines, respectively. The black lines show the result for equal mass planets. The other colors are for different mass ratios. Ejections become more common and collisions less frequent as the planetary mass ratio departs from unity.

#### Collisions

Collisions leave a single, larger planet in orbit around the star. The energy in the center-of-mass frame of the two planets is much smaller than the binding energy of a giant planet. Therefore, we model the collisions as completely inelastic and assume that the two giant planets simply merge together while conserving total momentum and mass. Under this assumption, we have calculated the distributions of orbital parameters for the collision products. The final orbit has a semi-major axis between the two initial semi-major axes, a small eccentricity, and a small inclination. While collisions between planets may affect the masses of extrasolar planets, a single collision between two massive planets does not cause significant orbital migration or eccentricity growth if the planets are initially on low-eccentricity, low-inclination orbits near the Hill stability limit.

#### Ejections

When one planet is ejected from the system, the less massive planet is nearly always ejected if $`\frac{m_<}{m_>}0.8`$, where $`m_<`$ and $`m_>`$ refer to the masses of the less and more massive planets, respectively. Since the escaping planet typically leaves the system with a very small (positive) energy, energy conservation sets the final semimajor axis of the remaining planet slightly less than

$$\frac{a_f}{a_1}\frac{a_2m_<}{a_1\left(m_1+m_2\right)+\left(a_2a_1\right)m_1}.$$
(1)
Thus, the ejection of one of two equal mass planets results in the most significant reduction in the semi-major axis, but is limited to $`\frac{a_f}{a_{1,\mathrm{init}}}0.5`$. The remaining planet acquires a significant eccentricity, but its inclination typically remains small. The eccentricity distribution for the remaining planet is not sensitive to the sum of the planet masses, but depends significantly on the mass ratio. In Fig. 2 (right) we show the cumulative distributions for the eccentricity after a collision for different mass ratios in different colors. While any one mass ratio results in a narrow range of eccentricities, a distribution of mass ratios would result in a broader distribution of final eccentricities. Also note that there is a maximum eccentricity which occurs for equal mass planets. Thus, the two planet scattering model predicts a maximum eccentricity of about 0.8 independent of the distribution of planet masses. While this compares favorably with the presently known planets, future observations will certainly test this prediction.

#### Stargrazers

In a small fraction of our numerical integrations ($`3\%`$) one planet underwent a close encounter with the central star. Due to the limitations of the numerical integrator used, the accuracy of our integrations for the subsequent evolution of these systems can not be guaranteed. Moreover, these systems could be affected by additional forces (e.g., tidal forces, interaction with the quadrupole moment of the star) that are not included in our simulations and would depend on the initial separation and the radius of the star. Nevertheless, our simulations suggest that for giant planets with initial semimajor axes of a few AU, the extremely close pericenter distances necessary for tidal circularization around a main sequence star (leading to the formation of a 51-Peg-type system) are possible, but rare.

### 2.3. Discussion

We now consider whether the the two planet scattering model could produce a distribution of eccentricities consistent with that of known extrasolar planets with periods longer than 10 days. In Fig. 3 (left) we show the cumulative distribution for the eccentricity of the remaining planet in our simulations. Here we have chosen a simple distribution of planet masses, $`P(m)m^1`$ for $`1M_{\mathrm{Jup}}m10M_{\mathrm{Jup}}`$, which is consistent with the mass distribution of the observed extrasolar planets (Tabachnik & Tremaine 2001). The red dashed line is just for systems which resulted in an ejection, while the green dotted line includes systems which resulted in a collision. Both are reasonably close to the observed distribution shown by the solid black line. It should be noted that the mass distribution used for this model has two cutoffs. While the $`10M_{\mathrm{Jup}}`$ cutoff is supported by observations, the lower cutoff is not constrained by present observations and is unlikely to be physical. Alternative mass distributions differing for $`m1M_{\mathrm{Jup}}`$ could alter this result. A detailed comparison to observations would require the careful consideration of observational selection effects, as well as the unknown initial mass distributions. Nevertheless, we find that the two-planet scattering model is able to reproduce the eccentricity distribution of the known extrasolar planets for plausible mass distributions.

## 3. Three Planets, Equal Masses

We have begun performing simulations of unstable systems containing three equal-mass planets using the conservative Bulirsch-Stoer integrator in Mercury (Chambers 1999). The results presented here are based on $`10^3`$ numerical integrations. For three-planet systems, we set the radius of each planet at the beginning of the integration. For determining when collisions occurred, each planet was assigned a radius of $`R_i=R_{\mathrm{J}up}`$. A planet was removed from the system if its separation from the central star became greater than $`1000`$AU (ejected) or less than $`r_{\mathrm{min}}=R_{}`$ (collided with star).

For three-planet systems, a single ejection or collision results in a system containing two planets. The remaining two planets may either be dynamically stable or undergo a subsequent ejection or collision to end with a single planet. Each run was terminated when two of the planets were removed by being ejected, colliding with another planet, or colliding with the star, or when the maximum integration time was reached $`t_{\mathrm{max}}=2\times 10^6`$ years.

Our numerical integrations were performed for a system containing three planets, each with a mass ratio $`m/M10^3`$. The initial semi-major axes are fixed at $`1`$, $`1.39`$, and $`1.92`$ AU, so that the ratio of orbital periods is $`P_3/P_2=P_2/P_1=1.63`$. Each planet is initially placed on a nearly circular orbit with an eccentricity of either $`0.0`$ or $`0.05`$ and an inclination of $`\pm 1^{}`$. All remaining angles (longitudes and phases) were randomly chosen between $`0`$ and $`2\pi `$.

The typical time for instability to develop is $`10^6`$ years. In the majority of cases, two planets were left in a stable configuration after one planet was ejected or two planets collided. In Fig. 3 (right) we show the distributions of the eccentricity of the inner planet after a collision or ejection of one planet. It is in reasonable agreement with the observed distribution for extrasolar planets. Both the eccentricity distribution (Fig. 3, right) and the corresponding distribution of the inclination of the inner planet (Fig. 4, left) are in qualitative agreement with those presented by Marzari & Weidenschilling (2002, based on slightly different initial conditions). In Fig. 4 (right) we show the cumulative distribution for $`\mathrm{\Delta }\varpi =\left|\left(\mathrm{\Omega }_2+\omega _2\right)\left(\mathrm{\Omega }_1+\omega _1\right)\right|`$, the angle between the periapses of the two planets remaining after one planet has been ejected from the system. We do not find the two remaining planets to have preferentially aligned periapses (cf. Malhotra 2002).

## 4. Conclusions

Dynamical interactions between planets of unequal mass reduces the frequency of collisions as compared to scattering between equal-mass planets. In particular, planet-planet scattering can reproduce the observed distribution of eccentric orbits using a plausible distribution of mass ratios. To determine if collisions are sufficiently rare will require a more careful comparison to the observations. Independent of the distribution of planet masses, planet-planet scattering by two planets initially on nearly circular orbits produces a maximum eccentricity of $`0.8`$. While this compares favorably with the presently known extrasolar planets, future detections will test this prediction of the two-planet scattering model. Alternatively, the observed properties of extrasolar planets could also be produced as a result of dynamical instabilities developing in three-planet systems.

##### Acknowledgments.

We are grateful to Eugene Chiang, Renu Malhotra, and Scott Tremaine for valuable discussions. E.B.F. acknowledges the support of the NSF graduate research fellowship program and Princeton University, and thanks the Theoretical Astrophysics group at Northwestern University for hospitality. F.A.R. and K.Y. acknowledge support from NSF grant AST-0206182. Our computations were supported by the National Computational Science Alliance under Grant AST980014N and utilized the SGI/Cray Origin2000 supercomputers at Boston University.

## References

Chambers, J.E. 1999, MNRAS, 304, 793

Chambers, J.E., Wetherill, G.W., & Boss, A.P. 1996, Icarus, 119, 261

Ford, E.B., Havlickova, M., & Rasio, F.A. 2001, Icarus, 150, 303

Gladman, B. 1993, Icarus, 106, 247

Lin, D.N.C. & Ida, S. 1997, ApJ, 447, 781

Malhotra, R. 2002, ApJ, L33

Marzari, F. & Weidenschilling, S.J. 2002, Icarus, 156, 570

Rasio, F.A. & Ford, E.B. 1996, Science, 274, 954

Rasio, F.A., Tout, C.A., Lubow, S.H., & Livio, M. 1996, ApJ, 470, 1187

Tabachnik, S. & Tremaine, S. 2001, astro-ph/0107482

Weidenschilling, S.J. & Marzari, F. 1996, Nature, 384, 619

Wisdom, J. & Holman, M. 1991, AJ, 102, 1528