text:
# Functional Correlation Approach to Operational Risk in Banking Organizations1footnote 11footnote 1The views presented in this paper are those of the authors and do not necessarily represent models or policies of Dresdner Bank AG.

## I Introduction

Risk management has become increasingly important in financial institutions over the last decade. Since the publication of JP Morgan’s RiskMetricsTM RM  in the nineties, Risk Management and Risk Control departments in banks have grown significantly in size and importance. The task is to fulfill regulatory requirements, to add transparency about a bank’s risk profile by a quantitative assessment of risks, to develop the necessary IT-solutions which allow to process the huge amount of data of a bank, and, finally, to integrate this information in a risk-return (RoRAC = Return on Risk-Adjusted Capital) based steering process of the bank. Ultimately, a proper risk management and risk control process is recognized by rating agencies and investors so that shareholder value is added to the bank.

Banks first focused on controlling potential losses due to market fluctuation, such as changes in the S&P 500 stock index, changes in interest and currency exchange rates, which is termed market risk. Internal market risk models are nowadays rather matured and accepted by regulators for the calculation of the required capital to be held as buffer against such losses. In contrast to these elaborated statistical models for market risks, credit risks (i.e., risks due to defaulted obligors) have to be covered by simply 8% capital of the bank’s risk-weighted assets. Implicitly, this charge also includes other risks such as operational risks. Since the New Basel Accord on Capital Adequacy issued by the Basel Committee on Banking Supervision in February and September 2001 BISIIa ; BISIIb ; BISIIc , known as Basel II, it is clear that regulators will demand banks to hold equity capital against operational risks explicitly.

A common industry definition of operational risk (OR) is the risk of direct or indirect losses resulting from inadequate or failed internal processes, people and systems or from external events Frachot01 . See Vandenbrink  for a practice-oriented introduction to the issue. Possible OR-risk categories are Vandenbrink : (i) human processing errors, e.g., mishandling of software applications, reports containing incomplete information, or payments made to incorrect parties without recovery, (ii) human decision errors, e.g., unnecessary rejection of a profitable trade or wrong trading strategy due to incomplete information, (iii) (software or hardware) system errors, e.g., data delivery or data import is not executed and the software system performs calculations and generates reports based on incomplete data, (iv) process design error, e.g., workflows with ambiguously defined process steps, (v) fraud and theft, e.g., unauthorized actions or credit card fraud, and (vi) external damages, e.g., fire or earthquake.

Thinking of theses categories as “operational risk processes” it is clear that there are functionally defined dependencies between individual processes, which all together bring a big organization to work. Consider the following example for illustration: a system error leads to an incomplete data import into a risk calculation engine, resulting in a wrong calculation of risk figures, and eventually to a human decision error by the trader, who closes a possibly profitable position unnecessarily to reduce a risk which in fact does not exist.

In the end misleading or lagging information, or system and workflow failures will always result in financial loss for a bank. Indeed, practitioners have recognized these dependencies in operational risk events and mandated units like the internal audit and risk control departments to control processes for the bank, and generated functions like a Chief Operating Officer (COO) to optimize them. Operational risk error trees between the above categories have been formalized in Vandenbrink  in more detail.

Since the mid nineties financial markets have also attracted physicists in academia. One of the main reason is that financial time series exhibit several statistical peculiarities, many of them being common to a wide variety of different markets and instruments. As such they could possibly be “universal”, i.e., independent of market details like instruments, country, and currency, and be the signature of collective phenomena in financial markets (see Bouch0 ; Bouch1 ; Bouch2 ; Sornette1 ; Stanly1  and references therein). Collective phenomena have been widely studied in physics in the context of phase transitions. Collective phenomena are often responsible for insensitivity of overall system behavior to details of an underlying dynamics. Specifically at phase transitions they give rise to power-law behavior, scale-invariance and self-similarity. Similar properties of financial time series might therefore well be understood as a consequence of agents in a market acting collectively.

Bringing together ideas from physics about collective phenomena and best industry practice for risk measurement the present paper details a possible statistical approach to determine the necessary equity capital to be held by banks to cover losses due to operational risks. In physical terms our model resembles a lattice gas with heterogeneous, functionally defined couplings. In such a description, bursts and avalanches of process failures correspond to droplet formation associated with a first order phase transition.

The paper is organized as follows. In Sect. II we describe the Value-at-Risk concept for risk management and control. In Sect. III we describe the approaches discussed in the context of Basel II for operational risk measurement. A new approach based on functionally dependent correlations giving rise to collective behavior is introduced in Sect. IV. Finally, Sect. V summarizes our results.

## II The Value-at-Risk Concept

Risk management in banks is based on diversification, hedging and equity capital as loss buffer. The bank charges its customers a premium for its risks so that (expected) losses in one market segment are on average compensated by profits in others. Other risks, especially market and increasingly credit risks, are hedged (insured) via the derivative market. Unexpected losses, which are not diversified or hedged, are covered by the bank’s equity capital. How much capital a bank needs to cover its risks is determined by the so-called “Value-at-Risk” (VaR). VaR can be defined as the worst loss in excess of the expected loss that can happen under normal market conditions over a specified horizon $`T`$ at a specified confidence level $`q`$. More formally, VaR measures the shortfall from the $`q`$-quantile of the loss distribution in excess of the expected loss, $`EL`$, within the time period $`T`$ discounted at the risk-free rate $`r`$ to time $`t=0`$:

$$\text{VaR}_{q,T}=\left(Q_q[L(T)]EL\right)\text{e}^{rT},$$
(1)
where the $`q`$-quantile worst case loss, $`Q_q[L(T)]`$, is defined at confidence level $`q`$ through

$$\text{Prob}\left(L(T)>Q_q[L(T)]\right)=1q.$$
(2)
As indicated in (1), VaR depends on the confidence level $`q`$ and the risk horizon $`T`$. The choice of these parameters depends on the application. If VaR is simply used to report or compare risks, these parameter can be arbitrarily chosen, as long they are consistent. If, however, VaR is used as a basis for setting the amount of equity capital, the parameters must be chosen with extreme care: the confidence level must reflect the default probability of the bank within the risk horizon, and the risk horizon must be related to the liquidation period of risky assets, recovery time of ill-functioning processes, or, alternatively, to the time period necessary to raise additional funds. This explains why regulators have chosen a high confidence level of 99% and a 10-day horizon to determine the minimum capital level for market risks. For credit risks and capital allocation, banks choose $`q`$ and $`T`$ even higher about 99.95% and one year, respectively.

In the financial industry there exist established statistical models for market and credit risk. Statistical models for operational risk start now to be discussed in the risk management community, especially in the context of Basel II BISIIb . Whereas internal market risk models are already recognized by regulators and are also used in banks for capital allocation, regulators are much more critical about internal statistical models for credit and operational risk. This is clearly less related to the mathematical complexity — although credit and operational risks are more difficult to model than market risks — but to problems with respect to input data which are much harder to validate than in the case of market risk.

## III Industry Standards for Operational Risk Measurement

The Basel Committee for Banking Supervision has proposed three alternative approaches to operational risk measurement BISIIc : The “Basic-Indicator Approach (BIA)”, the “Standardized Approach (SA)”, and the “Advanced Measurement Approach (AMA)”. In the BIA the required capital for operational risk is determined by multiplying a single financial indicator, which is gross income (interest, provision, trading, and other income) by a fixed percentage (called the $`\alpha `$-factor). The SA differs from the latter in that banks are allowed to choose business line specific weight factors, $`\beta _k`$, for the gross income indicator, $`I_k`$, of the $`k^{\mathrm{th}}`$ business line. The total regulatory capital charge, $`RC`$, is the simple sum of the capital required per business line,

$$RC=\underset{k}{}\beta _k\times I_k.$$
(3)
The weight factors $`\alpha `$ and $`\beta _k`$ are calibrated such that the required regulatory capital for operational risk would be 17 - 20% of the current regulatory capital on bank average standards.

The AMA consist of three sub-categories: The “Internal Measurement Approach (IMA)”, the “Loss Distribution Approach (LDA)” and the “Scorecard Approach (SCA)”. It is a more advanced approach as it allows banks to use external and internal loss data as well as internal expertise.

In the IMA the required capital is calculated as the sum over multiples of the expected loss per OR-risk category/business line cell

$$RC=\underset{i,k}{}\gamma _{ik}\times EL_{ik},$$
(4)
where $`i`$ is the risk category and $`k`$ the business line. The expected loss is quantified as the product of the annual OR-event probability, an exposure indicator per business line and risk category, and the loss percentage per exposure. All parameter estimates have to be disclosed to the supervisors. Since the $`\gamma `$-factor is computed on an industry based distribution, it will be possible to adjust the capital charge by a risk profile index, which accounts for the bank’s specific risk profile compared to industry.

The LDA and SCA are very similar as both approaches are based on a statistical VaR-model. Details of the LDA approach are outlined in ref. Frachot01 . In both approaches the bank estimates for each risk category/business line cell the probability distributions of the annual event frequency and the loss severity (= exposure $`\times `$ loss fraction per exposure). The difference between the LDA and the SCA is that in the former only internal or external historical loss data are used for estimating the distribution functions. In addition to this, banks are also allowed to apply expert knowledge to estimate the distribution functions in the SCA. This is a forward looking approach. It is particularly suited for operational risk, as processes that have failed are usually changed; hence historical loss data could provide potentially misleading information. Even if banks have an exhaustive internal database of losses, it can hardly be considered as representative of extreme losses. Hence, expert assessments and external loss databases are necessary. The problem with the latter data source is that external historical losses must be scaled to fit the balance sheet of the bank (it must be possible that the losses can occur in the bank).

Popular choices for the loss severity distribution functions are the lognormal, Gamma, Beta, Weibull distribution. Common choices for the loss frequency distribution function are the Poisson or negative binomial distribution. In a top-down approach different OR-risk categories are assumed to be independent. For each business unit and for each OR-category the OR-loss is simulated in a Monte-Carlo simulation by drawing a realization $`N_{ik}`$ from the loss frequency and sampling $`N_{ik}`$ realizations of the loss severity $`X_{ik}^m`$ ($`m=1,\mathrm{},N_{ik}`$). The loss in such a sample is

$$L_{ik}=\underset{m=1}{\overset{N_{ik}}{}}X_{ik}^m.$$
(5)
Drawing a histogram of outcomes of $`L_{ik}`$ provides the loss distribution function per risk category/business line cell. The Value-at-Risk is read off from the tail in excess to the expected loss as described in Sect. II. Due to the assumption of statistical independence, the loss distribution can also be calculated analytically as the convolution product of the loss frequency and the loss severity distribution. The required capital for the bank as a whole can either be calculated as the simple summation of the capital charges across each of the risk category/business line cell. This is the method given by the Basel Committee on Banking Supervision in the Internal Measurement Approach. Or, the MC-sampling can be extended beyond the risk category/business line cell by $`L=_{i,k}L_{ik}`$, which takes diversification between the risk category/business line cell into account.

A critical point which concerns all presently discussed approaches is the correlation between OR-losses. In this paper our focus will be how correlations and dependencies between OR-risk events can be integrated in the LDA/SCA.

## IV Functional Correlation Approach for Operational Risk

Since Markowitz’s centennial work on portfolio theory Markowitz52 ; Markowitz59 , diversification and dependencies between risk events are modeled by the covariance of stochastic processes. Because empirically only the mean and the covariance of these processes are reliably determined from market data, it is common practice to choose correlated Gaussian white noise for modeling correlations. As a consequence the loss distribution is unimodal with frequent small losses and a few extreme losses, which — dependent on the distribution of the loss severity — are responsible for fat tails in the loss distribution. Collective losses or even crashes such as burst and avalanches of losses are not contained in this description.

A main point in this paper is that this stochastic dependency of risk events is not sufficient for all risk categories: one frequently also observes direct, functional and non-stochastic dependencies. Functional dependency between risk events is most pronounced in operational risk events. Processes in a (large) organization are usually organized so as to mutually support each other. Thus, if a process fails, this will usually be detrimental to other processes relying on receiving input or support of some sort from the failing process in question, so that they run a higher risk of failing as well. It therefore seems inadequate to model operational risk events individually per risk category/business line cell and aggregate losses afterwards over some covariance matrix, which would be the choice when approaching operational risks analogously to market risks. In the following we extend the LDA/SCA by taking the functional dependencies between processes into account.

We consider a simple two state model here, i.e., a processes can be either up and running or down. For the process corresponding to the OR-event $`i`$ we designate these states as $`n_i=0`$ and $`n_i=1`$, respectively. In following we will skip the business line index $`k`$ for simplicity.

The interest is in obtaining reliable estimates of the statistics of processes that are down at any time and of the statistics of losses incurred at any time. As the loss severity incurred by a given process going into the down state may vary randomly from event to event, solving the latter problem requires convolving the statistics of down-events with the loss severity distribution related to the process failures.

The reliability of individual processes will vary (randomly) across the set of processes, and so will the degree of functional interdependence. These random heterogeneities constitute an element of quenched disorder, whereas the loss severities incurred by down processes constitute an element of annealed disorder as they are (randomly) determined anew from their distribution each time a process goes down. An appealing feature for the modeler of operational risks therefore is the independence of the dynamic model of the interacting processes and the loss severity model (i.e. the estimate of the PDFs of loss severity incurred by individual process failures.) A typical assumption for the latter is to take them as being distributed according to a log-normal distribution with suitable parameters for means and variances, which we will choose in the following.

### IV.1 Dynamics

To motivate the dynamics of the functional approach, note that all processes need a certain amount of “fueling” or support in order to maintain a functioning state for the time increment $`tt+\mathrm{\Delta }t`$ within the risk horizon, $`t[0,T)`$ (think of human resources, information, input from other processes, etc.). Here, only the generic features of the model shall be outlined. Hence, the increment $`\mathrm{\Delta }t`$ is chosen such that all processes can fully recover within this time interval, i.e., the state $`n_i`$ of each process can flip each time step. For practical applications in banks, one would model the recovery process more carefully: specific death-period after the failure of the $`i^{\mathrm{th}}`$ process would be considered, and one would differentiate between process failures being discovered and adjusted up to a certain cut-off time, e.g., end-of-day, at which a process would have been completed Vandenbrink . These features are not generic and can only be discussed related a specific OR-event under consideration.

We denote by $`h_i(t)`$ the total support received by process $`i`$ at time $`t`$, and choose it to take the form

$$h_i(t)=\vartheta _i\underset{j}{}w_{ij}n_j(t)+\eta _i(t).$$
(6)
That is, it is composed of (i) the average total support $`\vartheta _i`$ that would be provided by a fully operational network of processes (in which $`n_i(t)=0`$ for all $`i`$). This quantity is (ii) diminished by support that is missing because of failing processes which normally feed into the process in question; (iii) lastly, there are fluctuations about the average which we take to be correlated Gaussian white noise with — by proper renormalizing $`\vartheta _i`$ and $`w_{ij}`$ — zero mean and unit variance. Correlated Gaussian noise is introduced to model equal-time cross correlations between OR-risk categories in analogy to the approach proposed by the Basel Committee for Banking Supervision for credit risk,

$$\eta _i(t)=\sqrt{\rho }Y(t)+\sqrt{1\rho }ϵ_i(t),$$
(7)
where $`Y(t)𝒩(0,1)`$ is a common factor for all OR-risk categories with equal-time correlation coefficient $`\rho `$, and the $`ϵ_i(t)𝒩(0,1)`$ are idiosyncratic terms.

Note that non-linear effects could be included by modifying (6) to $`h_i(t)=\vartheta _i_jw_{ij}n_j(t)_{j,k}w_{ijk}n_j(t)n_k(t)\mathrm{}+\eta _i(t)`$. Note also that, as in credit risk modeling, the common factor $`\sqrt{\rho }Y(t)`$ could further be decomposed into sector-contributions $`\sqrt{\rho }Y(t)_k\beta _{ik}Y_k(t)`$ so as to describe more complicated equal-time correlations. To keep this treatment transparent, we will present the formalism without these extensions.

Process $`i`$ will fail in the next time instant $`t+\mathrm{\Delta }t`$, if the total support for it falls below a critical threshold. By properly renormalizing $`\vartheta _i`$, we can choose this threshold to be zero, thus ($`\mathrm{\Theta }`$ is the step-function: $`\mathrm{\Theta }(x)=1`$ for $`x0`$ and 0 else)

$`n_i(t+\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{\Theta }(\vartheta _i+{\displaystyle \underset{j}{}}w_{ij}n_j(t)`$ (8)
$`\sqrt{\rho }Y(t)\sqrt{1\rho }ϵ_i(t)).`$
The losses incurred by process $`i`$ are then updated according to

$$L_i(t+\mathrm{\Delta }t)=L_i(t)+n_i(t+\mathrm{\Delta }t)X_{t+\mathrm{\Delta }t}^i,$$
(9)
where $`X_{t+\mathrm{\Delta }t}^i`$ is randomly sampled from the loss severity distribution for process $`i`$. Note that the process dynamics is independent of assumptions concerning their loss severity distributions within the present model.

One can integrate over the distribution of idiosyncratic noises to obtain the conditional probability for failure of process $`i`$ given a configuration $`n(t)=\{n_i(t)\}`$ of down-processes and a realization of the common factor $`Y(t)`$ at time $`t`$,

$`n_i(t+\mathrm{\Delta }t)_{n(t),Y(t)}`$
$`\text{Prob}\left(n_i(t+\mathrm{\Delta }t)=1|n(t),Y(t)\right)`$
$`=\mathrm{\Phi }\left({\displaystyle \frac{\mathrm{\Phi }^1(p_i)+_jw_{ij}n_j(t)\sqrt{\rho }Y(t)}{\sqrt{1\rho }}}\right).`$
Here $`\mathrm{\Phi }(x)`$ denotes the cumulative normal distribution. Note that we have set $`\vartheta _i=\mathrm{\Phi }^1(p_i)`$, where $`p_{i,\mathrm{\Delta }t}p_i`$ is the unconditional expected probability for process failures within the time-increment $`\mathrm{\Delta }t`$. This is consistently justified by setting $`n_j(t)=0`$ for all $`j`$ and $`\rho =0`$ in Eq. (IV.1). Note that up to the functional term “$`_jw_{ij}n_j(t)`$” this approach corresponds to the approach adapted by the Basel Committee for Banking Supervision for credit risk BISIIc .

The couplings $`w_{ij}`$ can be determined by considering the transition probabilities, $`p_{ij,\mathrm{\Delta }t}p_{ij}`$, for process $`i`$ failure within the time-increments $`\mathrm{\Delta }t`$, given that in the configuration at time $`t`$ process $`j`$ is down, whereas all other processes are running, and $`Y(t)=0`$. Introducing the shorthand $`𝒞_j`$ for this configuration, we can write

$`p_{ij}`$ $`=`$ $`\text{Prob}\left(n_i(t+\mathrm{\Delta }t)=1|𝒞_j\right)`$ (11)
$`=`$ $`\mathrm{\Phi }\left({\displaystyle \frac{\mathrm{\Phi }^1(p_i)+w_{ij}}{\sqrt{1\rho }}}\right).`$
This leads to

$$w_{ij}=\sqrt{1\rho }\mathrm{\Phi }^1(p_{ij})\mathrm{\Phi }^1(p_i).$$
(12)
Analogous identities would be available for determining higher order connections $`w_{ijk}`$, if nonlinear effects were taken into account. Note that the probabilities for process failure depend only on the increment $`\mathrm{\Delta }t`$ and not on the time $`t`$ due to the stationarity of the dynamics.

To illustrate how these parameter are fixed in practice, consider the following. Either from a historical loss database, or from an expert assessment the following two questions must be answered per OR-risk category and business line:

1. What is the expected period, $`\tau _i`$, until process $`i`$ fails for the first time in a fully operative environment, and
2. given that only process $`j`$ has failed, what is the expected period, $`\tau _{ij}`$, for process $`i`$ to fail also?

Noting that with $`\text{Prob}(\mathrm{failure}\mathrm{at}z\mathrm{\Delta }t)=\left(1p_i\right)^{z1}p_i`$ one finds that

$$\tau _i=\underset{z=1}{\overset{\mathrm{}}{}}z\mathrm{\Delta }t\left(1p_i\right)^{z1}p_i=\frac{\mathrm{\Delta }t}{p_i},$$
(13)
and analogously,

$$\tau _{ij}=\underset{z=1}{\overset{\mathrm{}}{}}z\mathrm{\Delta }t\left(1p_{ij}\right)^{z1}p_{ij}=\frac{\mathrm{\Delta }t}{p_{ij}}.$$
(14)
These identities express the $`p_i`$ and $`p_{ij}`$ in terms of estimated average times of failure, and are used to fix the model parameters completely. Note that according to (11) $`p_{ij}`$ can be interpreted as a non-equal time correlation for process failures.

Note also that, incidentally, the dynamics (8) resembles that of a lattice gas (defined on a graph rather than on a lattice), the $`n_i`$ denoting occupancy of a vertex, the $`w_{ij}`$ interactions, and $`\vartheta _i`$ taking the role of chemical potentials regulating a-priori occupancy of individual vertices. The present system is heterogeneous in that (i) the $`\vartheta _i`$ vary from site to site, (ii) the couplings $`w_{ij}`$ have a functional rather than regular geometric dependence on the indices $`i`$ and $`j`$ designating the vertices of the graph. Moreover, in the physics context, one usually assumes noise sources other than Gaussian so that cumulative probabilities are described by Fermi-functions rather than cumulative normal distributions as above. The quantitative difference is minute, however.

The model dynamics as such cannot be solved analytically for a general heterogeneous network. We shall resort to Monte-Carlo simulations to study its salient properties. The main qualitative features can, however, also be observed in the simplified situation of a homogeneous network consisting of identical processes, having the same connectivity at each node. A mean-field analysis of such a simplified situation will be given as well. As the presence of the common factor expressed by the $`\rho `$-term in Eq. (7) would influence only quantitative details of the system’s behavior, we will further present the analysis without correlation to the common factor by setting $`\rho 0`$.

### IV.2 Key Features

Key features of the collective behavior of networks of interacting processes can easily be anticipated either directly from a discussion of the dynamic rules, or from the analogy with the physics of lattice gasses.

Consider a network in which the unconditional probabilities for process failures, $`p_i`$, are small, but process interdependence is large and consequently conditional probabilities for process failures, $`p_{ij}`$, are sizeable. In such a situation, spontaneous failure of individual processes may induce subsequent failures of other processes with sufficiently high probability so as to trigger a breakdown of the whole network. If, on the other hand, process interdependence remains below a critical threshold value, individual spontaneous failures will not have such drastic consequences, and the whole network will remain in a stable overall functioning state.

Of particular interest for the risk manager is the case, in which process interdependence is low enough to make a self-generated break-down of the network extremely unlikely, but parameters are nevertheless such that a stable overall functioning state of the network coexists with a phase in which nearly the complete network is in the down state (two phase coexistence). In such a situation, it may be external strain which can induce a transition from a stably functional situation to overall breakdown. Analogous mechanisms are believed to be responsible for occasional catastrophic breakdowns in bistable ecosystems Nature01 .

With increasing unconditional probabilities for process failures it becomes meaningless to distinguish between an overall functioning, and a non-functioning phase of the network, two–phase coexistence ceases to exist — as in (lattice) gasses — at a critical point.

### IV.3 Simulations

In the following we validate some of our intuitions about global network behavior using Monte-Carlo simulations.

The Monte-Carlo dynamics can either be conceived as parallel dynamics (all $`n_i`$ are at each time step simultaneously updated according to (8) or (IV.1)), or as (random) sequential dynamics (only a single $`n_i`$ is (randomly) selected for update according to (8) in any given time-step, in which case the time increment must scale with the number $`N`$ of processes in the net as $`\mathrm{\Delta }tN^1`$).

For the analysis of operational risks, losses are accumulated during a Monte-Carlo simulation of the process dynamics over the risk horizon, $`T`$. Runs over many risk horizons then allow to measure loss distribution functions for individual processes within the network of interacting processes, or of business units or the full network by appropriate summations.

For the simulations, we choose a random setting, i.e., unconditional failure probabilities are taken to be homogeneously distributed in the interval $`[0,p^{\mathrm{max}}]`$ and we determine random conditional failure probabilities as $`p_{ij}=p_i(1+\epsilon _{ij})`$, with $`\epsilon _{ij}`$ homogeneously distributed in $`[0,\epsilon ^{\mathrm{max}}]`$, which fixes the ratio $`(p_{ij}/p_i)^{\mathrm{max}}`$.

Fig. 1 shows a situation where a functional network coexists with a situation in which the network is completely down, and parameters are such that spontaneous transitions between the phases are not observed during a simulation. The upper track shows the loss record of a system initialized in the “all down” state whereas the lower track exhibits the loss record of the same network initialized in the fully functioning state.

The loss distribution for the functional network is unimodal with a bulk of small losses and a fat tail of extreme losses, which are driven by the loss severity distribution.

By increasing the functional interdependence at unaltered unconditional failure probabilities, the functioning state of the network becomes unstable. A spontaneous transition into the ‘down’ state is observed during a single run of 50000 Monte-Carlo steps. Two interesting features about this transition to complete breakdown deserve mention: (i) the time to breakdown can vary within very wide limits (we have not attempted to measure the distribution of times to breakdown and its evolution with system parameters such as range of values for conditional failure probabilities), (ii) there are no detectable precursors to the transition; it occurs due to large spontaneous fluctuations carrying the system over a barrier, in analogy to droplet formation associated with first order phase transitions.

We should like to emphasize that realistically the system dynamics after an overall break-down of a process network would no longer be the spontaneous internal network dynamics: recovery efforts would be started, increasing support for each process by a sufficient amount such as to reinitialize the network in working order.

Repeated spontaneous transitions in both directions can be observed only in a rather small network (Fig. 2). The corresponding loss distribution will be bimodal.

Fig. 3 illustrates the principle of a strain-simulation. In each case, the system, if in the operational state, is repeatedly put under external strain by turning off 5 randomly selected functioning processes every 1000th time step, and letting the system evolve under its internal dynamics thereafter. Such a disturbance can either trigger a breakdown of the system or not. In the former case, if the system is found fully down 1000 time steps later, it is reinitialized in the fully operational state and once more disturbed 1000 steps later.

One observes that the operational low-loss phase which we have seen in Fig. 1 to coexist with the non-operational phase is resilient against disturbances of the kind described above. The same is true if $`(p_{ij}/p_i)^{\mathrm{max}}`$ is increased to 2.7. At $`(p_{ij}/p_i)^{\mathrm{max}}=2.8`$ an external strain succeeds once during the simulation to trigger breakdown of the net, whereas at $`(p_{ij}/p_i)^{\mathrm{max}}=2.9`$ breakdown under external strain of the given strength is the regular response of the system (with a few exceptions and occasional spontaneous recoveries).

Of interest to the risk manager in the end are the total losses accumulated over a risk horizon, $`T`$,

$$L(T)=\underset{i}{}L_i(T),$$
(15)
more specifically, the corresponding probability density function. Fig. 4 presents such distribution of accumulated losses for a network that remains operational throughout the simulation. For this simulation we have chosen $`T=365\mathrm{\Delta }t`$. The loss distribution has an extended tail (barely visible on the scale of the data, with a 99.5% quantile at 1400, and the largest aggregated loss observed during the simulation over a time span of $`T`$ at 3400, i.e. by a factor of more than 3 larger than the expected loss for the chosen risk horizon $`T`$.

### IV.4 Mean-field Solution

The principle dynamic properties of the model exist independently of its heterogeneous nature. To exhibit systematic relations between dynamic features and model parameters, we shall elucidate them in the simplified setting of a homogeneous network of identical processes, which we analyze within a mean-field approximation.

For this we assume a homogeneous coupling, $`w_{ij}w_0/z`$ and $`p_{ij}p_w`$, where $`z`$ is the coordination number of the graph (taken to be identical at each vertex), and replace time–dependent quantities in Eq. (IV.1) by long-time stationary averages, $`n_j(t)n_j(t)`$. In a homogeneous system, averages are independent of the process index $`i`$ and time, such that $`p_ip`$ and $`n_j(t)n`$. This gives the mean-field equation (without correlation to the common factor, $`\rho 0`$)

$`n=\mathrm{\Phi }\left(\mathrm{\Phi }^1(p)+w_0n\right).`$ (16)
Depending on $`p`$ and the average coupling strength, $`w_0`$, this equation has either one unique solution or three solutions, with one unstable solution at intermediate $`n`$, and two stable solutions, $`n`$ 0 or 1. Figure 5 shows stable and unstable solutions as functions of $`w_0`$ and $`p_w`$ for a given value of $`p`$.

The phase diagram (Fig. 6) summarizes regions in the $`p`$$`w_0`$ plane, where operational and non-operational phases of the network coexist showing limits of stability of the low-loss and high-loss solution. For $`p`$ exceeding a critical value $`p_c0.0218`$ there is a unique disordered phase with relatively large values of $`n`$.

For a sufficiently small unconditional failure probability $`p`$ an initially running process network will for weak functional dependence remain in the running state, despite of spontaneous individual process failures. Such a network is in a functioning state.

For stronger functional correlation the functional state of the network becomes unstable. Fluctuations in the number of processes that are down at any given time can trigger a burst or avalanche of failures — a collective phenomenon corresponding to droplet formation associated with a first order phase transition.

For intermediate degrees of functional dependence, the network allows for two (meta)stable states, an operational and a non-operational network. While spontaneous fluctuations in the number of failing processes will in large networks fail to trigger transitions between the two states, external strain may well do, as demonstrated in the strain simulations above. Seen from the functioning side, the closer the parameters are to an instability line, the smaller the strain needed to trigger avalanches of failures leading to the fully defunct state of the system. This behavior is quantified by the susceptibility $`\chi `$ (Fig. 7). With support for each process decreasing by an amount $`\delta h<0`$, the fraction of down processes will change by

$$\delta n\chi (w_0,p)\delta h$$
(17)
with

$$\chi (w_0,p)=\frac{1}{1w_0\mathrm{\Phi }^{}(\mathrm{\Phi }^1(p)+w_0n))},$$
(18)
where $`\mathrm{\Phi }^{}(x)=\frac{1}{\sqrt{2\pi }}\mathrm{exp}[x^2/2]`$. Note that the susceptibility is proportional to the sensitivity of the fraction of down-processes with respect to the unconditional probability of process failure.

Triggering complete failures becomes easier close to instability lines for two reasons: (i) the susceptibility diverges as instability lines are approached; (ii) the unstable solution $`n_u`$ is closer to the stable equilibrium solution; external strain only needs to push the system beyond $`n_u`$ to destabilize the functioning state.

Such a regime could therefore be dangerous for a network of mutually supporting processes in a bank. Due to the long periods in one of the metastable states, the bank would not necessarily realize the potential of big losses due to bursts and avalanches of process failures. With a basically unchanged process setup the network could collapse and cause significant losses, either due to external strain or rare fluctuations of internal dynamics. Owing to the stability of the metastable states, the bank will then have to spend a lot of efforts in order to bring the network back to a functional state, which will cause additional costs.

## V Conclusion

In this paper we have outlined how ideas from physics of collective phenomena and phase transitions can naturally be applied to model building for operational risk in financial institutions. Our main point was that functional correlations between mutually supportive processes give rise to non-trivial temporal correlation, which could eventually lead to the collective occurrence of risk event in form of burst, avalanches and crashes. For risks associated to process failure (operational risks) a functional dependence seems to be the appropriate way for modeling sequential correlations.

From the physics point of view, the appropriate model is rather simple, being a heterogeneous variant of the well studied lattice gas model. Despite the heterogeneities, it has a first order phase transition (driven by average interaction strength) at sufficiently low average a-priori probability of process failures, showing coexistence between an overall functioning state (gas) and a state of catastrophic breakdown (liquid). As the a-priori probability of process failures is increased the first order transition ends at a (liquid/gas) critical point.

One of the most critical lessons for Risk Control from our analysis is the possible metastability of networks of interacting processes: The bank would not necessarily realize the potential of big losses due to bursts and avalanches of process failures, as there are no detectable precursors to such transitions. With a basically unchanged process setup the network could collapse and cause significant losses, either due to external strain or rare fluctuations of internal dynamics. Owing to the stability of the metastable states, the bank will then have to spend a lot of efforts in order to bring the network back to a functional state, which will cause additional costs. To assess the metastability banks have to perform stress tests.

It should be noted that realistically the system dynamics after an overall break-down of a process network would no longer be the spontaneous internal network dynamics: recovery efforts would be started, increasing support for each process by a sufficient amount such as to reinitialize the network in working order.

In a forthcoming publication we will also show that the random walk model for financial time series commonly used in banks can naturally be extended to incorporate functional dependencies leading to collective effects. This will lead to models which, while bearing some resemblance with agent-based models of markets (see, e.g., agents ) are different from them in other respects tobepub .

Acknowledgments We would like to thank U. Anders and G. J. van den Brink of Dresdner Bank for valuable discussions and practical input.


format: text

text:
# Deconstructing Dense Coding

\[

## Abstract

The remarkable transmission of two bits of information via a single qubit entangled with another at the destination, is presented as an expansion of the unremarkable classical circuit that transmits the bits with two direct qubit-qubit couplings between source and destination.

\]

Quantum dense coding enables Alice to communicate two bits of classical information by sending Bob a single physical qubit, which is maximally entangled with another qubit already in his possession. She does the trick by applying one of four unitary transformations to her member of the entangled pair, thereby converting the state of the pair into one of four mutually orthogonal two-qubit states. Bob can learn which state it is after receiving the second member of the pair.

What is surprising is that Alice appears to act on only a single qubit, thereby providing Bob with two bits of information by sending him only one appropriately prepared qubit. But this way of telling the tale downplays a second interaction that takes place before the curtain even rises on the official story. That earlier interaction is required to create the entanglement between the qubits that Alice and Bob initially share.

The full story remains surprising even with this added prologue, but the real surprise is that the entangling interaction, essential for the transmission of the two bits, can take place before Alice has even chosen the bits she wishes to communicate to Bob. What the story really demonstrates is the remarkable ability of entangled states to store interaction in a highly fungible form that need not be cashed in until the need arises.

I have made a similar point about quantum teleportation, showing explicitly how the missing interaction that makes the difference between a routine classical circuit and a quantum “miracle”, is buried in the interaction that produces a crucial shared entangled pair before the state to be teleported need even have been formed. Because teleportation and dense coding both exploit preexisting shared entanglement to facilitate communication with surprisingly little additional interaction, one might expect there to be a similar circuit-theoretic deconstruction of dense-coding. But because there is no direct mapping from one protocol to the other — teleportation involves three qubits and dense coding four — it is not obvious from the expansion in  of classical state swapping into quantum teleportation, how dense coding might arise from an expansion of the classical circuit that communicates two bits of information by means of two direct qubit-to-qubit interactions.

In this note I show how to do this. The construction is given in Fig. 1. The generalization from qubits to $`d`$-state systems is given in Fig. 2 and Eqs. (1)-(3).

If the initial state of the four qubits in Fig. 1(a) is $`|x|y|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0}`$ (reading from top to bottom on the left) then the effect of the two cX (cNOT) gates is to transform it into $`|x|y|x|y`$. This automates a classical procedure by which Alice, who possesses the upper two qubits, can communicate two classical bits of information to Bob, who possesses the lower two.

To go from this undramatic classically transparent procedure to quantum dense coding we first expand the cX on the right into quantum components, beginning with the fact (Fig. 1(b)) that $`\text{X}=\text{H}\text{Z}\text{H}.`$ This is useful because we wish to eliminate, or at least disguise, the direct coupling on the left between Alice and Bob’s lower qubit. Because the operator Z is diagonal in the computational basis, it is immaterial whether Z acts on a control qubit immediately before or immediately after a cX. So since cX is its own inverse we can expand Fig. 1(b) to Fig. 1(c), and then move the paired cX and Hadamard gates to the extreme left and right, as shown in Fig. 1(d). The goal of eliminating the direct coupling between Alice and Bob’s lower qubit can now be achieved by noting that the two cX gates on the left of Fig. 1(d) are equivalent to the three cX gates on the left of Fig. 1(e), since both sets, acting on the computational basis, leave the control qubits unaltered, while applying X to the lowest qubit if and only if the states of the two control qubits differ. But since Bob’s qubits both start on the left in the state $`|\mathrm{\hspace{0.17em}0}`$, and X acts as the identity on $`\text{H}|\mathrm{\hspace{0.17em}0}=\frac{1}{\sqrt{2}}\left(|\mathrm{\hspace{0.17em}0}+|\mathrm{\hspace{0.17em}1}\right)`$, the leftmost cX in Fig. 1(e) always acts as the identity and can be dropped from the circuit.

The result, Fig. 1(f), is an automated dense coding circuit. The two gates on the left convert $`|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0}`$ into the maximally entangled state $`\frac{1}{\sqrt{2}}\left(|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0}+|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}1}\right).`$ The upper member of the entangled pair is then acted on by $`\text{X},\text{Z},\text{Z}\text{X}`$ or no transformation at all, depending on whether the state of the upper two qubits is $`|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}1}`$, $`|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}0}`$, $`|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}1}`$, or $`|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0}`$. The two gates on the extreme right then transform the resulting entangled state of the two lower qubits (one of the four states of the “Bell basis”) back to whichever computational basis state of the upper two qubits gave rise to it.

A generalization of the dense-coding protocol from qubits to $`d`$-state systems has recently been given by Liu et al. In the corresponding generalization of the circuit-theoretic derivation the cX operator becomes the controlled bit rotation,

$$\text{cX}:|x|y|x|yx,0x,y<d,$$
(1)
where $``$ denotes addition modulo $`d`$, the Hadamard transformation H becomes the quantum Fourier transform

$$\text{H}:|y\frac{1}{\sqrt{d}}\underset{0z<d}{}e^{2\pi izy/d}|z,$$
(2)
and the controlled-$`Z`$ operation becomes

$$\text{cZ}:|x|ye^{2\pi ixy/d}|x|y.$$
(3)
One easily verifies that

$$(\text{H})_2(\text{cX})_{12}=(\text{cZ})_{12}^{}(\text{H})_2$$
(4)
and therefore

$$\text{cX}_{12}=(\text{H})_2^{}(\text{cZ})_{12}^{}(\text{H})_2.$$
(5)
Fig. 2 extends the identities of Fig. 1 to $`d`$-state systems. The only difference in the diagrams is that the unitary gates are no longer self-inverse, and must be distinguished from their adjoints. Fig. 2(a) shows two direct couplings by controlled bit rotations (1) that take $`|x|y|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0}`$ into $`|x|y|x|y`$, $`0x,y<d`$. Fig. 2(b) introduces the identity (5). A controlled bit rotation and its compensating inverse are introduced in Fig. 2(c). The replacement of the two cX gates on the left of Fig. 2(d) by the two cX and one $`\text{cX}^{}`$ gates on the left of Fig. 2(e) is clearly valid for controlled bit rotations, and the cX gate on the left of Fig. 2(e) can be dropped since $`\text{H}|\mathrm{\hspace{0.17em}0}`$ is invariant under arbitrary bit rotations.

Fig. 2(f) is the $`d`$-state version of dense coding. The two gates on the left produce the entangled state

$$\frac{1}{\sqrt{d}}\underset{0z<d}{}|z|z.$$
(6)
The two gates in the middle transform (6) by the action (or inaction) of the $`\text{cX}^{}`$ and $`\text{cZ}^{}`$ gates on the member of the entangled pair in Alice’s possession. The two gates on the right act on the pair after both members are in Bob’s possession, transforming its state into that product of Alice’s two computational-basis states that governed the two controlled operations in the middle.

These circuit-theoretic deconstructions of dense coding (and the corresponding deconstructions of teleportation in ) back into elementary classical circuits, illustrate the role of entanglement as interaction-in-advance-of-need, by explicitly tracing its origin back to a direct classical interaction. They have the pedagogical virtue of requiring no algebraic scratchwork whatever (except for the confirmation of (5) for $`d`$-state systems) to verify that the quantum circuits act as advertised.

This work was supported by the National Science Foundation, Grant PHY0098429.


text:
# Conductance fluctuations at the quantum Hall plateau transition

## Abstract

We analyze the conductance fluctuations observed in the quantum Hall regime for a bulk two-dimensional electron system in a Corbino geometry. We find that characteristics like the power spectral density and the temperature dependence agree well with simple expectations for universal conductance fluctuations in metals, while the observed amplitude is reduced. In addition, the dephasing length $`L_\mathrm{\Phi }T^{1/2}`$, which governs the temperature dependence of the fluctuations, is surprisingly different from the scaling length $`L_{sc}T^1`$ governing the width of the quantum Hall plateau transition.

Conductance fluctuations are probably one of the most prominent features occurring in the mesoscopic world. For metallic (diffusive) systems at low enough temperature these so-called universal conductance fluctuations (UCF) are characterized by an amplitude $`e^2/h`$, independent on sample size and material, and, they are generally believed to be well understood Lee et al. (1987).

Similar conductance fluctuations are also observed in mesoscopic quantum Hall (QH) systems, however, the underlying physics still remains an open question. Interpretations of different experiments cover modifications of UCF in high fields G.Timp et al. (1987); Geim et al. (1992); Hohls et al. (1999), tunnelling between opposite edge states through bulk inhomogeneities Simmons et al. (1991); Main et al. (1994); Bykov et al. (1996), influences of charging effects Cobden and Kogan (1996); Cobden et al. (1999), and networks of compressible-incompressible regions Machida et al. (2001). Common to all the previous experiments is a sample geometry in the form of a (Hall-)bar, the importance of edge channels, and a typical sample width $`w`$ of 1-3 $`\mu `$m. Transport in edge channels and the coupling of both edges due to single impurities is strongly relevant, making the applicability of a general UCF-theory problematic.

Here we will address the conductance fluctuations in the quantum Hall regime due to transport through the disordered bulk of a two-dimensional electron system (2DES). In order to avoid any edge effects we use a Corbino geometry. The large width $`w=6\mu `$m, considerably above the elastic scattering length $`l_{\mathrm{el}}<1\mu `$m, suppresses transport through individual impurities. Analyzing the temperature dependence of small fluctuations superimposed on the conductance peak of the QH plateau transition we find a behavior very common to UCF in metallic systems, however, with a considerably reduced absolute value of the UCF amplitude.

The samples used for this work are based on modulation doped GaAs/AlGaAs heterostructures. Additional impurities (Be or Si) added into the 2DES Ploog (1987); Haug et al. (1987) yield an enhanced short range scattering and a small elastic scattering length. Throughout this paper we present data measured on a sample with an electron mobility $`\mu =2`$ m2/Vs, an electron density $`n=2.110^{15}`$ m-2 and a density of Be-impurities $`n_{Be}=210^{14}`$ m-2. We confirmed our results with measurements on other samples, containing both Si- and Be-impurities, with mobilities ranging up to $`\mu =12`$ m2/Vs. Using annealed AuGe-Ni contacts the samples were patterned into Corbino geometry with a ring width $`w=6\mu `$m, and an inner radius $`r=60\mu `$m.

The samples were mounted on the cold finger of a dilution refrigerator with a base temperature $`T<20`$ mK and positioned into the center of a superconducting solenoid. We measured the current with a Lock-In technique for fixed AC voltage (9 Hz) as a function of magnetic field and temperature. We carefully checked for heating effects by varying the amplitude of the excitation voltage at the same temperature. Excitation amplitudes adapted to the different temperature regimes were used to avoid heating at the lowest temperatures and still ensure high enough resolution of the fluctuations at elevated temperatures. The voltage was fixed to $`1\mu `$V for $`T<50`$ mK, $`2\mu `$V for $`50T120`$ mK, and $`5\mu `$V for $`T>120`$ mK.

Due to their low mobilities our samples show broad quantum Hall plateaus with vanishing conductance $`G=I/V`$ around integer filling factors $`\nu =nh/eB`$. Near half integer filling factors the conductance is nonzero with a peak value of the order of $`10e^2/h`$, corresponding to diagonal conductivities $`\sigma _{xx}=(1/2\pi )\mathrm{ln}(1+w/r)G0.15e^2/h`$, with an actual value between 0.1 and 0.3$`e^2/h`$ (depending on sample and filling factor). Here we concentrate on the plateau transition $`\nu =21`$ which is broadest in magnetic field. The conductance at this transition is shown in Fig. 1a for different temperatures. It is dominated by the well understood conductance peak which arises when the Fermi energy sweeps through the delocalized electronic states near the center of a spin split Landau band.

Superimposed onto this conductance peak are small reproducible conductance fluctuations, see Fig. 1b. For further analysis they first have to be separated from the peak form $`\overline{G(B)}`$ obtained by a polynomial smoothing. Only features on scales $`\mathrm{\Delta }B>50`$ mT are included into $`\overline{G(B)}`$. The final fluctuations $`\delta G(B)=G(B)\overline{G(B)}`$ extracted from the peak are shown in Fig. 1c. Distinct reproducible features can be clearly observed, they become smoother with rising temperature and their amplitude shrinks.

Considering, that the electronic transport in the QH plateau transition with extended states at the Fermi energy may be regarded as quasi metallic, it is worthwhile to compare the observed conductance fluctuations with UCFs in mesoscopic 2D metals, which are governed by the minimum of either the dephasing length, $`L_\varphi =\sqrt{D\tau _\varphi }`$, or the thermal length, $`L_T=\sqrt{\mathrm{}D/k_BT}`$ Lee et al. (1987). Here $`D`$ is the diffusion coefficient and $`\tau _\varphi `$ is the phase coherence time.

For the low temperatures considered here ($`T<1`$ K) it was shown that in high magnetic fields electron-electron scattering dominates $`\tau _\varphi `$ in a 2DES realized in a semiconductor Brandes (1995). For a typical conductivity $`\sigma _c0.5e^2/h`$ near the critical point of the QH plateau transition it is estimated to $`\tau _\varphi \mathrm{}/(k_BT)`$ Brandes et al. (1994); Polyakov and Samokhin (1998). In this framework the dephasing length $`L_\varphi `$ approximately equals the thermal length $`L_T`$ and a distinction between $`L_\varphi `$ and $`L_T`$ becomes unnecessary.

The dephasing length $`L_\varphi (T)`$ influences the temperature dependence of both the amplitude of the conductance fluctuations, $`\sqrt{\delta G^2}`$, and their correlations. This fact can be formalized when using a temperature dependent correlation function averaged over the magnetic field, $`F(\mathrm{\Delta }B)=\delta G(B)\delta G(B+\mathrm{\Delta }B)_B`$. Defining a correlation field, $`B_c(T)\mathrm{\Phi }_0/L_\varphi (T)^2`$, the correlation function then embraces two limits Lee et al. (1987): It is constant for $`\mathrm{\Delta }BB_c`$, i.e. $`F(\mathrm{\Delta }B)=F(0)=\delta G^2`$, reduced to $`F(B_c)=F(0)/2`$ for $`\mathrm{\Delta }B=B_c`$, and it follows $`F(\mathrm{\Delta }B)1/\mathrm{\Delta }B^2`$ for $`\mathrm{\Delta }BB_c`$.

UCFs can be analyzed in terms of the power spectral density (PSD), defined as the Fourier transform of the correlation function $`F(\mathrm{\Delta }B)`$

$$𝒫(f_B)=\frac{1}{2\pi }F(\mathrm{\Delta }B)\mathrm{exp}(i\mathrm{\hspace{0.17em}2}\pi f_B\mathrm{\Delta }B)𝑑\mathrm{\Delta }B$$
(1)
with $`f_B`$ the magnetic frequency. It has been shown Schäfer et al. (1999) that using the PSD yields a more accurate analysis of the data in an experimental situation like ours, where we are only able to measure fluctuations on a magnetic field scale $`\mathrm{\Delta }B<0.05`$T, and, where we are additionally limited to a finite range of the magnetic field ($`B=6.47.1`$ T).

In Fig. 2a the PSD $`𝒫(f_B)`$ of our experimentally measured conductance fluctuations, $`\delta G(B)=G(B)\overline{G(B)}`$, is shown for several temperatures. Each point within the graph is the result of averaging $`\mathrm{ln}(𝒫(f_B))`$ over a frequency interval of $`\mathrm{\Delta }f_B=10`$ T-1. The data follow an exponential decay indicated by the straight lines. The observed saturation of the $`T=160`$ mK data for frequencies $`f_B60`$T-1 is caused by measurement noise which is constant with frequency and therefore always dominates the PSD at high $`f_B`$. For $`T>200`$ mK the frequency interval, where the PSD does not sink into the measurement noise, becomes too small for a trustworthy analysis.

The exponential decay of the PSD observed in Fig. 2a can be fitted by a simple model

$$𝒫(f_B,B_c(T))=𝒫_0e^{2\pi B_c(T)f_B}.$$
(2)
Here $`𝒫_0`$ is a temperature independent prefactor. The corresponding correlation function $`F(\mathrm{\Delta }B)`$ fulfills all predictions for the limits of $`F(\mathrm{\Delta }B)`$. Additionally, the temperature dependence of $`\delta G^2=F(0)`$ extracted from this PSD follows the expected behavior for a 2D metal, $`\delta G^2L_\varphi ^2(T)1/B_c(T)`$ Lee et al. (1987). A least square fit shown by the straight lines in Fig. 2a now enables us to determine the correlation field $`B_c(T)`$ plotted in Fig. 2b. For $`T50`$ mK the correlation field follows a linear dependence $`B_cT`$ as expected for $`B_c=\mathrm{\Phi }_0/L_\varphi ^2`$ and $`L_\varphi T^{1/2}`$. The corresponding temperature dependent dephasing length $`L_\varphi `$ is shown on the right axis of Fig. 2b. Down to a temperature $`T=50`$ mK it is still well below the width of the sample, $`w=6\mu `$m. The observed saturation of $`L_\varphi `$ at temperatures below $`T=50`$ mK is presumably caused by a decoupling of the electron temperature from the bath temperature due to microwave heating.

Using the experimentally measured PSD we can now directly calculate the amplitude square $`\delta G^2`$ of the UCFs. However, since our procedure of extracting $`\delta G`$ from the peak disregards fluctuations on a scale $`\mathrm{\Delta }B>50`$mT we have to restrict ourselves to frequencies $`f_B`$ above a cutoff frequency $`f_g>(50\mathrm{mT})^1`$. This defines a reduced fluctuation amplitude square

$$\delta G^2_{\mathrm{\Delta }B<(1/f_g)}=2_{f_g}^{\mathrm{}}𝑑f_B𝒫(f_B)$$
(3)
where mainly fluctuations on small magnetic field scales $`\mathrm{\Delta }B<1/f_g`$ are taken into account. For a large enough $`f_g`$ effects of the smooth background substraction are negligible and we can safely extract an experimental UCF amplitude to be compared to theoretical calculations.

Fig. 3 shows the measured reduced amplitude square of the fluctuations as defined in Eq. 3 for $`f_g=40`$T-1. The data are shown for the complete magnetic field interval ($`6.4\mathrm{T}<B<7.1`$T) and for small $`B`$-intervals on both flanks of the peak around the points of half peak height ($`6.4\mathrm{T}<B<6.6`$T and $`6.9\mathrm{T}<B<7.1`$T). The experimental data are compared with the expectations from UCF theory, which using Eqs. (2) and (3) yields $`\delta G^2_{\mathrm{\Delta }B<1/f_g}=𝒫_0/(\pi B_c)\mathrm{exp}(2\pi B_cf_g)`$ with a critical field $`B_c`$ linearly depending on $`T`$. As shown in Fig. 3 the temperature dependence of all three data sets can be described by UCF theory.

Using this reduced amplitude as given in Eq. 3 or equivalent the experimentally determined PSD as given in Eq. 2 we can extrapolate from the raw experimental data to the full fluctuation amplitude $`\sqrt{\delta G^2(T)}0.02(e^2/h)(1/\sqrt{T(\mathrm{K})})0.1(e^2/h)L_\varphi (\mu `$m). This value can be compared to the fluctuation amplitude

$$\sqrt{\delta G^2(T)}\frac{e^2}{h}\left(\frac{2\pi r}{w}\right)^{1/2}\frac{L_\varphi (T)}{w}$$
(4)
for a ring shaped metallic sample with radius $`r`$ and width $`w`$ Lee et al. (1987). An extension of UCF theory to the presence of Landau quantization does not predict any change in the universal fluctuation amplitude Xiong and Stone (1992); Maslov and Loss (1993); Khmelnitskii and Yosefin (1994); Xiong et al. (1997). With the parameters $`w=6\mu `$m and $`r=60\mu `$m for our sample we expect $`\sqrt{\delta G^2}(e^2/h)L_\varphi (\mu `$m), one order of magnitude larger than the experimental value deduced above. This discrepancy mirrors that the 2DES at the QH plateau transition can be described to some extend like a 2D metal but is not a metal. In fact, the conductivity $`\sigma _{xx}<e^2/h`$ does not allow the perturbative treatment to first order of $`(e^2/h)/\sigma _{xx}`$ used in UCF theory Lee et al. (1987); Xiong and Stone (1992) for metallic systems. However, numerical studies of the two-point conductance $`G`$ of a rectangular sample at the transition from the Hall insulator $`G=0`$ to the lowest Hall plateau $`G=1`$ with a Chalker and Coddington network model Cho and Fisher (1997); Jovanovic and Wang (1998) also observe a reduced amplitude $`\sqrt{\delta G^2(T)}0.3e^2/h`$ compared to $`\sqrt{\delta G^2(T)}e^2/h`$ for high filling.

Finally, the temperature dependence of the fluctuations amplitude can be used to compare the length scale $`L_\varphi L_TT^{1/2}`$ governing the fluctuations with the length $`L_{sc}(T)T^{1/z}`$ appearing in the scaling behavior of the QH plateau transition, $`z`$ is the so-called dynamical scaling exponent (for a review see Huckestein (1995)). For temperatures where the sample width $`w`$ exceeds $`L_{sc}`$ the transition width of the QH plateau transition scales as $`\mathrm{\Delta }BL_{sc}^{1/\gamma }T^\kappa `$ with $`\kappa =1/z\gamma `$. For lower temperatures $`\mathrm{\Delta }B`$ was experimentally found to saturate Koch et al. (1991). This saturation could be attributed to finite size scaling occuring for $`w<L_{sc}`$ which allowed a direct measurement of the localization length in the QH plateau transition.

Usually $`L_{sc}`$ is identified as the dephasing length $`L_\varphi `$. However, it was noted that for $`L_\varphi >w`$ one would expect large mesoscopic fluctuations Liu and Das Sarma (1994). In contrast, no fluctuations Koch et al. (1991) or only small fluctuations Koch et al. (1992) were observed experimentally at the lowest temperatures. As a consequence, the observed saturation of $`\mathrm{\Delta }B`$ was claimed to be caused by other effects such as external heating of the 2DES Liu and Das Sarma (1994). In our experiments we can use the temperature dependence of the fluctuation amplitude as a thermometer for the 2DES. As can be seen from the temperature dependence of the correlation field in Fig. 2b and the fluctuation amplitude $`\delta G^2`$ in Fig. 3 an electron temperature $`T_e`$ well below 50 mK is achieved. For the lower conductance at the flanks, where unwanted microwave heating is reduced, $`\delta G^2`$ even follows the model down 30 mK (Figs. 3 and 4). In contrast, as shown in Fig. 4 the width $`\mathrm{\Delta }B`$ of the conductance peak only follows a scaling law $`\mathrm{\Delta }BT^\kappa `$

for $`T150`$ mK and saturates for $`T<100`$ mK. Interpreted as a saturation decoupling of electron temperature from the temperature in the mixing chamber one would find a minimum electron temperature $`T_e=70`$ mK inconsistent with the observed temperature dependence of the fluctuations. This demonstrates that the observed transition width saturation is indeed an effect of the sample size. Additionally, using the generally accepted value $`\gamma =2.3`$ of the critical scaling exponent and applying the relation $`\kappa =1/z\gamma `$ for the critical exponents of the QH plateau transition one finds $`L_{sc}1/T`$ Huckestein (1995). This deviates substantially from the temperature dependence of $`L_\varphi 1/\sqrt{T}`$. Both observations demonstrate that size scaling and conductance fluctuations are governed by distinct length scales implying that $`L_{sc}`$ depends on different mechanisms than $`L_\varphi `$, an experimental finding in agreement with recent theoretical predictions Polyakov and Samokhin (1998).

In conclusion, we have measured the conductance fluctuations of a 2DES in the quantum Hall regime. Their dependence on temperature and magnetic field can be described by standard UCF-theory for metals, however, their amplitude is considerably reduced. In addition we have shown that temperature dependence of the quantum Hall plateau transition and that of the fluctuations are governed by different length scales.

The samples used for our experiments were grown by K. Ploog at the Max Planck Institut für Festkörperforschung. We thank F. Evers, F. Kuchar, D. Polyakov, and L. Schweitzer for stimulating discussions.


text:
# Stability and eccentricity of periodic orbits for two planets in a 1:1 resonance.

## 1 Introduction

In many extra-solar planetary systems discovered recently, the observed Doppler shift of the emitted light is well described by assuming that the central star is moving in a Keplerian elliptic orbit due to its gravitational interaction with a single Jupiter-size planet. If residuals are present in a least squares fit to the data after possible chromospheric fluctuations in the star have been taken into account, these residuals signal the presence of an additional planet, or possibly several planets (Fischer et al, 2000) (Marcy et al 2001a, ). There is, however, an important exception when such a star also travels on a Keplerian orbit even though there are two planets orbiting it. We have in mind Lagrange’s celebrated solution of the three-body problem (Lagrange, 1772), for which he won the prize of the Royal Academy of Science of Paris. In this solution, each body moves on a Keplerian orbit with a common plane, period, eccentricity and focus that is located at their center of mass, in such a manner that at all times the relative positions of these bodies form the vertices of an equilateral triangle of variable size, (see Fig. 1). This solution will be discussed in section 2. It is therefore interesting to consider the possible occurrence of such a 1:1 resonance in extra-solar planetary systemm with two Jupiter-size planets, particularly in view of the recent discovery of a remarkable 2:1 resonance of two such large planets in GJ876 (Marcy et. al. 2001b, ). To be relevant to current astronomical discoveries, however, it is necessary that Lagrange’s solution be stable in the range of observed masses and eccentricities. In the past, linear stability analyses have been carried out that were primarily focused on the restricted three-body problem (where one of the three masses vanishes), which is applicable to the study of the motion of asteroids, but these results can also be extended to the general Lagrange solution as will be discussed in section 3. In this paper we present the results of a numerical analysis of the nonlinear stability domain of Lagrange’s solution as a function of the masses of the three bodies and the eccentricity of the common elliptical orbits (see section 4). In practice there will be deviations from Lagrange’s solution, and the interesting question arises as to whether there are distinct characteristics which would distinguish these deviations from other types of perturbations due to additional planets. While variations in the eccentricity and major axis of the approximate elliptical orbit are common to all perturbations caused by the presence of a second planet, one of the most distinguishing feature of a slightly off-resonance but stable Lagrange solution is that the the rotation rate of the axis of the elliptic orbit of the star is much smaller than for other types of perturbations, (see Fig. 2).

For example, for a $`1\%`$ deviation in the position of the lighter planet, assumed to be about a Jupiter mass, we find that it takes about $`800`$ periods for a complete revolution of the major axis, (see Fig. 2). In contrast, a fit to the recently discovered 2:1 resonance in GJ876 indicates that the major axis of the planets should complete a revolution in about $`53`$ periods of the heavier planet (Nauenberg, 2002).

Lagrange’s solution is not the only 1:1 resonance for the three-body problem. We also consider another solution that we found to be stable in the domain of masses and eccentricities relevant to extra-solar planetary systems (Laughlin 2002, ). In this case the two lighter bodies (planets), and the heaviest body (star) have different orbits, and therefore the Doppler shift data should be readily distinguishable from the case of a single planet. The characteristic feature of this solution is that at each half period the three bodies are aligned, and if the heavier planet is in a nearly circular orbit, the lighter planet moves in a highly eccentric elliptical orbit, (see Fig. 3). A periodic alignment of the two planets and the star is also a characteristic of other resonance solutions, as in the case of the 2:1 resonance in GJ876(Laughlin and Chambers, 2001) (Rivera and Lissauer, 2001) (Lee and Peale, 2002) (Nauenberg, 2002). The evolution of this configuration when the planets are slightly off resonance is shown in Fig. 4, which illustrates the clockwise rotation of the major axis of the eccentric orbit.

## 2  Lagrange’s solution of a three-body problem

In Lagrange’s solution (Lagrange, 1772) of a three-body problem, each body travels on a separate elliptic orbit with a common period, eccentricity and focus which is located at their center of mass, in such a way that these bodies are always at the vertices of an equilateral triangle of variable size. An example is illustrated in Fig.1, where the masses are in the ratio $`1.0:0.7:0.2`$, and the common eccentricity is $`ϵ=.5`$, showing the equilateral triangle for the relative positions of the three bodies at apogee. Apart from the overall orientation of the system, these parameters uniquely describe the solution. In this solution the relation between variables is somewhat different from the corresponding ones in the two-body problem. For example, the common frequency $`\omega `$ or period $`P=2\pi /\omega `$ of the motion is given by the relation

$$\omega =\sqrt{Gm(1+ϵ)^3/R^3},$$
(1)
where $`G`$ is Newton’s constant, $`m=m_1+m_2+m_3`$ is the total mass, $`m_1`$ is the mass of the heaviest body (star), $`m_2`$, $`m_3`$ are the masses of the lighter bodies (planets), $`ϵ`$ is the common eccentricity, and $`R`$ is the maximum size of the equilateral triangle on which the three bodies are located. Then the major axis $`a_i`$ of each of the elliptic orbits is

$$a_i=\frac{\sqrt{m_j^2+m_k^2+m_jm_k}}{m(1+ϵ)}R,$$
(2)
where the subscripts $`i,j`$ and $`k`$ are permutations of the the integers $`1,2,3`$, while the mean of the velocity of the star at the maximum and minimum distance from the foci of the ellipse is

$$K=(2\pi G/P)^{1/3}\frac{\sqrt{m_2^2+m_3^2+m_2m_3}}{m^{2/3}\sqrt{1ϵ^2}}.$$
(3)
If the planetary masses $`m_2`$ and $`m_3`$ are small compared to the mass of the star, $`m1`$, then $`a_2a_3a`$, where $`a=R/(1+ϵ)`$ and $`\omega \sqrt{Gm/a^3}`$ as in the corresponding two-body problem. In principle, any fit to the data with a single planet of mass $`m_p`$ can also be attributed to two planets which, according to Eq. 3 for $`K`$, have masses satisfying the relation $`m_p=\sqrt{m_2^2+m_3^2+m_2m_3}`$, which is somewhat less than the sum of the masses of the two planets. In practice, however, it is unlikely that extra-solar planetary systems would occur in an exact 1:1 resonance, and therefore the presence of a second planet manifests itself in the occurrence of residuals in a single Keplerian orbit fit to the data.

## 3 Linear stability analysis

The first linear stability analysis of Lagrange’s solution was carried out by Routh for the special case of circular orbits (Routh, 1875) 111In this paper, Routh incorrectly attributed Lagrange’s solution to Laplace.. Assuming that the attractive forces between the bodies depends on the relative distance $`r`$ as $`1/r^\kappa `$, Routh demonstrated that Lagrange’s solution was stable provided that the masses satisfied the inequality

$$\gamma <\frac{1}{3}(\frac{3\kappa }{1+\kappa })^2,$$
(4)
where

$$\gamma =\frac{(m_1m_2+m_1m_3+m_2m_3)}{(m_1+m_2+m_3)^2}.$$
(5)
For gravitational interactions where $`\kappa =2`$, the constant on the right hand side 222 Routh found that he had been anticipated in this important result by someone called M. Gascheau in a thesis on mechanics, apparently left unpublished  of Eq. 4 is $`1/27`$, which implies that the masses of the two lighter bodies (planets), $`m_2`$ and $`m_3`$, must be much smaller than the mass of the heaviest body (star) $`m_1`$. For example, setting $`m_2=0`$, which correspond to the restricted three-body problem and applies to the motion of asteroids such as the Trojans, this inequality implies that $`m_3/(m_1+m_3)<.03852..`$ which has become known as Routh’s critical point. Neglecting quadratic terms in the mass ratios $`m_2/m_1`$ and $`m_3/m_1`$, Routh’s inequality becomes approximately $`(m_2+m_3)/m_1<1/27`$, indicating that the stability depends to a very good approximation only on the sum of the masses of the lighter bodies.

For the stable configurations, Routh obtained the frequencies $`\omega _1,\omega _2`$ and $`\omega _3`$ of the normal modes in the plane of the orbits, which for $`\kappa =2`$ are given by

$$\omega _1=\omega \sqrt{\frac{1}{2}(1+\sqrt{127\gamma )}},$$
(6)
$$\omega _2=\omega \sqrt{\frac{1}{2}(1\sqrt{127\gamma )}},$$
(7)
$$\omega _3=\omega ,$$
(8)
where $`\omega `$ is the fundamental Kepler frequency, Eq. 1, and he determined the corresponding amplitudes for these modes. While Routh did not discussed the stability with respect to deviations perpendicular to the plane of the orbit, it is straightforward to show that for this case the Lagrange orbits are stable. Remarkably , he also briefly considered second order deviations, which he remarked could “ultimately disturb the stability”, but this part of his analysis was incomplete, calling attention only to two of the possible commensurability or resonance relations, $`\omega _3=2\omega _2`$ and $`\omega _1=2\omega _2`$, for which $`\gamma =1/36`$ and $`\gamma =16/675`$, respectively.

For elliptic orbits, a linear stability analysis of Lagrange’s solution was not carried out until some 90 years later, when (Danby, 1964) numerically integrated the Floquet equations for the first order deviations from the Lagrange solution of the restricted three - body problem (see also (Bennett, 1965). Subsequently, a majority of stability studies have been confined to this special case, but as we shall see, the results can also be applied to the general solution (Marchal 1990), ). For a modern discussion of the stability of the Lagrange solution see (Siegel and Moser, 1971)

## 4 Nonlinear stability domain of Laplace’s solution

The nonlinear stability domain of Lagrange’s solution for the three-body problem that is presented here was obtained by integrating the equations of motion numerically, and determining whether, for a small initial deviation from the solution, the orbits were either confined or unconfined after a large number $`n`$ of periods. We found that when $`n`$ was increased from 400 to 800, there were no significant changes in the results except near the critical points discussed below, where we increased $`n`$ until no further changes occurred. For the deviations we consider small displacements of the position of one of the lighter bodies planets in the plane of the orbit and also perpendicular to this plane. As we shall see, except in the case when $`\omega _3=2\omega _2`$, only initial deviations in the plane of the orbit gave evidence for some of the expected nonlinear instabilities near the commensurability relations for the three fundamental frequencies, Eqs. 6-8.

Starting with a small deviation in the velocity of one of the lighter bodies of order $`1\%`$ in the direction perpendicular to the plane of the orbits, our numerical result for the nonlinear stability domain is shown in Fig. 5, where the unstable regions are indicated by small squares. In agreement with the linear stability analysis ( see section 3), the nonlinear stability domain depends only on the common eccentricity $`ϵ`$ of the orbits (horizontal axis), and the Routh parameter $`\gamma `$, Eq. 5 (vertical axis). This parameter is determined by the ratios $`m_2/m_1`$ and $`m_3/m_1`$ of the masses of the lighter bodies (planets) to the mass $`m_1`$ of the heaviest body (star) (vertical axis). We have verified this result by computing this domain for different fixed values for the ratio of the masses of the planets, without observing any changes when plotting the results with Routh’s parameter $`\gamma `$. For the particular case shown in Fig. 5, the computation was carried out for $`m_2=m_3`$, while the curves shown are the linear stability computations of Danby and Bennett (Danby, 1964) (Bennett, 1965), which were originally obtained for the restricted three-body problem where either $`m_2`$ or $`m_3`$ is set equal to zero.

This nonlinear stability domain looks surprisingly similar to the stability domain of the Mathieu equation (Grimshaw, 1990) with an additional nonlinear restoring or damping term. Indeed, for the restricted three-body problem, a linear stability analysis of the equations of motion, to first order in the common eccentricity $`ϵ`$, in a frame of reference rotating with the frequency of the orbits, yields a bifurcation at $`\omega _3=2\omega _2`$ which can be viewed as a parametric resonance between the fundamental orbital frequency for elliptic motion and the oscillation frequency of the first order deviations of the massless body. In the general case, we have seen that this bifurcation occurs for $`ϵ=0`$ at $`\gamma =1/36=.02777\mathrm{}`$. For the restricted problem, $`m_2=0`$, this corresponds to $`\mathrm{\Delta }_0=2\sqrt{2}/3.9428`$, where $`\mathrm{\Delta }=(m_1m_3)/(m_1+m_3)`$ and $`m_3/m_1.02944`$. We evaluated the instability domain to first order in the eccentricity $`ϵ`$, and found that it lies inside the wedge

$$\mathrm{\Delta }_0\delta ϵ\mathrm{\Delta }\mathrm{\Delta }_0+\delta ϵ,$$
(9)
where $`\delta =3/20\sqrt{2}.10606`$. This domain is shown in Fig. 6, together with the corresponding nonlinear domain where the orbit actually becomes unstable. These linear boundaries mark the onset of a bifurcation in which the orbits first become aperiodic filling out a confined region of space.

An enlargement of the stability domain near the critical point at $`ϵ=0`$ and $`\gamma =1/27.037037..`$ is shown in Fig. 7, which displays our results for the case that $`m_2=m_3`$. This critical point is often discussed in connection with the stability of the Lagrangian points for the restricted three body problem in a frame of reference rotating with the frequency of the circular orbits, and is applied to the study of the Trojan and other asteroids (Murray & Dermott, 1999). As in the previous case, however, this critical point marks only a bifurcation to an aperiodic but confined motion, while numerically we find that the nonlinear instability begins at a somewhat higher value of $`\gamma .0391`$, corresponding for $`m_2=0`$ to $`m_3/m_1.0425`$. A curve quadratic in the eccentricity fits well the boundary of the upper part of the nonlinear instability domain, as shown in Fig. 7, while a linear curve fits the lower part of the domain with a slope which is somewhat higher than the one which we calculated analytically in the linear approximation, Eq. 9 but in good agreement with the numerical linear stability results of (Danby, 1964) and (Bennett, 1965).

Up to now, we have considered only initial deviations in the direction perpendicular to the plane of Lagrange’s orbits. In Fig. 8 we show the nonlinear stability domain obtained by taking an initial a displacement $`dx=dy=.001`$ in this plane for one of the two lighter bodies. This result was obtained by fixing the ratio $`m_2/m_3=1`$, but the same results are obtained for other values of the ratios of these masses when the results are plotted as a function of the Routh parameter $`\gamma `$, Eq. 5. We see evidence in this figure for the instabilities at the nonlinear resonances $`\omega _1=2\omega _2`$ and $`\omega _3=3\omega _2`$, corresponding to $`\gamma =16/675.02370..`$ and $`\gamma =32/2187.01463..`$ Also, the domain of stability shrinks in the upper wedge with indications of additional resonances in this region. This is confirmed by an enlargement of this region, which shows the tail of the instability due to the resonance at $`3\omega _1=4\omega _2`$ corresponding to $`\gamma =576/16875.03413\mathrm{}`$ as well as other resonances (see Fig. 9).

## 5  Another 1:1 resonance solution

Lagrange’s solution is not the only stable 1:1 resonance system for the three-body problem. Another solution is illustrated in Fig. 3, in which the two lighter bodies or planets are traveling on two different orbits which are approximately elliptic in such a way that the central star and the planets are aligned when located at the maximum or minimum distance from the center of mass. We found these periodic orbits (Laughlin 2002, ) by a new method based on an expansion of the coordinates in a Fourier series with a common period, and we determined the Fourier coefficients of this expansion by finding minima of the action integral with respect to these coefficients by an iterative process (Nauenberg, 2001). A characteristic feature of these orbits is that when the heavier of the two planets is in a nearly circular orbit, the lighter planet is in a elliptical orbit with a very large eccentricity which rises slowly with increasing mass of the heavier planet, as shown in Fig. 10 .

In contrast to the Lagrange’s solution discussed previously, a small deviation from exact 1:1 resonance in this case leads to a relatively rapid rotation of the major axis of the elliptical orbit, as illustrated in Fig 3. This leads to characteristic modulations in the Doppler shift oscillations of the light emitted by the star. As in the case of Lagrange’s solution discussed in the previous section, we found that these orbits are stable over the range of masses relevant to extra-solar planetary systems, see Fig. 11.

## 6 Conclusions

The nonlinear stability domain for Lagrange’s solution of the three body problem shown in Fig. 8 indicates that there is a wide range of Jupiter-size planetary masses (including brown dwarfs) and eccentricities for which such solutions can exist in extra-solar planetary systems. For example, for an eccentricity of $`ϵ.6`$ the ratio of the total mass of the two planets to the mass of the star for which the solutions are stable is $`.004`$, except for a small region where nonlinear resonances occur. This mass correspond to $`4.2`$ Jupiter-masses for a one solar-mass star, while for smaller eccentricities, $`ϵ.2`$, there is a wedge of stable solutions for higher mass ratios up to approximately $`.04`$. In principle, any Doppler shift data that can be fitted under the assumption that only a single planet is orbiting the central star, can equally well be attributed to two planets orbiting the star according to Lagrange’s solution. In practice, however, it is very unlikely that two planets are in an exact 1:1 resonance, and therefore one expects to find residuals in the data which signal the present of a second planet. One of the main effects due to the perturbations caused by a second planet is the secular rotation of the major axes of the approximate Keplerian ellipses that characterize the orbits of the planets and the star. This effect, however, would not occur in the case of an exact 1:1 Lagrangian resonance, and it is strongly suppressed in the case that the resonance is approximate. It may be thought that in view of the greater number of degrees of freedom present with Lagrange’s solution, a better least squares fit to the data should be readily available. This is not the case, because Lagrange’s equations also allow for unphysical solutions where the mass of one or even both of the planets can have negative values, provided the total mass is positive. Indeed, in a preliminary attempt to obtain a least squares fit to data which show residuals, the optimal mass of the smaller of two Lagrange planets turned out to be negative, which ruled out this solution.

The stability domain for the type of 1:1 resonance solution discussed in section 5 is shown in Fig. 11, demonstrating that this solution also encompasses the possibility of two Jupiter-size planets orbiting a solar-size star. In this analysis we restricted the heavier planet to be in a nearly circular orbit, and found that the lighter planet is in a highly eccentric orbit with $`ϵ.8`$ , (see Fig. 4). If this restriction is relaxed, we also find similar stable solutions, and for equal masses the two planets exchange eccentricity when the major axis rotates through $`180`$ degrees.

In summary, it is likely that extra-solar planetary systems that have several Jupiter-size planets which are close enough to give rise to significant gravitational perturbations will be in resonance, because numerical investigations have shown that such systems can be stable. In such cases, the planets and the central star are periodically aligned. An interesting exception is the 1:1 resonance solution of Lagrange, where the planets and the star are located at all times on the vertices of an equilateral triangle of varying size. It would be very exciting if this solution, discovered by Lagrange 230 years ago, and realized thus far only in the motion of the Trojan and other asteroids in our solar system (Murray & Dermott, 1999), would be also present in the orbits of planets in extra-solar systems. Likewise a search should be undertaken to find also two planets in extra-solar systems which are in a 1:1 resonance of the type discussed in section 5, which does not occur in our system.

I would like to thank Richard Montgomery and Carles Simo for useful references to the vast literature on the Lagrange solution of a three-body problem. michael@mike.ucsc.edu.