text:
# AF embeddings and the numerical computation of spectra in irrational rotation algebras

## 1. Introduction

In this paper we address the problem of finding numerical approximations to the spectrum of a bounded linear operator on a complex, separable Hilbert space. This is a very difficult problem in general and there are various ways to attack it as well as various quantities that one may wish to use to approximate the spectrum (e.g. pseudospectra or numerical ranges). A very natural approach is to start with a given operator $`T`$ acting on a Hilbert space $`H`$ and “compress” it to a finite dimensional subspace of $`H`$. This compression is a finite dimensional matrix whose eigenvalues can hopefully be numerically computed and one further hopes that these eigenvalues will somehow approximate the spectrum of the infinite dimensional operator $`T`$. This approach has been studied by numerous authors (see, for example, ,  and the references therein) and a number of interesting results have been obtained for large classes of operators.

Though the present work grew naturally out of the author’s own study of the method described above (cf. \[4, Section 6\]) it turns out that an absurdly abstract approach to numerical approximations of spectra actually proves useful in at least one important case. The main result of this paper simply observes how some very specialized and technical work in C-algebra theory can be used to get excellent approximations to spectra of some important operators. In particular, this applies to discretized one-dimensional Schrödinger operators with quasiperiodic potential. See  for a nice treatment of the spectral theory of Almost Mathieu operators (up to 2001) and  (and its references) for a recent contribution to the more general case of quasiperiodic Schrödinger operators.

When presented abstractly the basic idea of this note becomes quite simple, but will require a bit of C-algebra theory to explain. The first basic fact we need is that all finite dimensional C-algebras have a special form – they are just finite direct sums of finite dimensional matrix algebras. We will also need the definition of an AF (“approximately finite dimensional”) algebra; A C-algebra $`C`$ is AF if there exist finite dimensional C-subalgebras $`C_1C_2\mathrm{}C`$ whose union is dense in $`C`$. The observation below is the trivial, but key, idea of this paper.

Observation: If $`TB(H)`$ is an element in an AF C-subalgebra of $`B(H)`$ then there exist finite dimensional matrices whose spectral quantities (e.g. pseudospectra or numerical ranges) approximate those of $`T`$.

Indeed, if $`TCB(H)`$ with $`C`$ an AF algebra then we can find operators $`T_n`$ such that $`TT_n0`$ (hence spectral quantities of $`T_n`$ are close to those of $`T`$) and each $`T_n`$ is contained in a finite dimensional C-algebra (hence can be identified with a finite dimensional matrix). There is one subtlety that must be mentioned. Namely, the actual spectrum of $`T_n`$ need not be close to the actual spectrum of $`T`$ (unless $`T`$ happens to be normal) but singular values, numerical ranges and pseudospectra will be close. Moreover, the actual spectrum of $`T`$ is equal to the intersection of all its pseudospectra and hence we still get reasonable approximations to the spectrum of $`T`$ by looking at pseudospectra of the $`T_n`$’s. In the case both $`T`$ and $`T_n`$ are normal (e.g. self-adjoint Schrödinger operators) then the actual spectra are close in a very strong sense; the Hausdorff distance between their spectra is bounded by $`TT_n`$.

AF-embeddability (i.e. studying which operators belong to AF algebras) has been considered by numerous authors but most of this work is too general to be of much help in finding explicit matrices which approximate given operators. For example, it is an easy consequence of the spectral theorem that every normal operator is contained in an AF-subalgebra of $`B(H)`$ but this gives no clue how to actually find the right matrix approximations. However, returning to the seminal paper on AF-embeddability (cf. ) one finds that not only are explicit matrix models provided (for operators in irrational rotation algebras) but computable and sharp rates of convergence can be proved.

In section 2 we will review the necessary aspects of the theory of continued fractions. We also state a few simple lemmas that will be needed later.

In section 3 we review the construction of Pimsner-Voiculescu which provides explicit matrix models for AF-embeddings of the irrational rotation algebras. We will not reproduce the proofs of any estimates as these can be found either in the original paper  or in \[5, Chapter VI\].

In section 4 we state and prove the main result for spectral approximations of operators in irrational rotation algebras. The proof really amounts to working out a few estimates as the hard part was accomplished more than 20 years ago in .

Finally, in section 5 we show how results of Haagerup-Rørdam can be used to give general “one-sided” rates of convergence. Though the results of this section are not as good as previous sections they have the advantage of always being numerically implementable.

## 2. Preliminaries

In this section we will review some basic facts about continued fractions and point out a few technical (but very elementary) inequalities that we will need later. Though there are a number of excellent books on continued fractions, a classic reference is . The facts we will need, however, can be found in any book on the subject. We also state the definition of pseudospectrum at the very end of this section.

Given an irrational number $`\theta (0,1)`$ there exist unique positive integers $`a_1,a_2,\mathrm{}`$ such that

$$\theta =\underset{n\mathrm{}}{lim}[a_1,\mathrm{},a_n]=\underset{n\mathrm{}}{lim}\frac{1}{a_1+\frac{1}{a_2+\frac{1}{\mathrm{}+\frac{1}{a_n}}}}.$$

Note that if $`\theta `$ is given then one can compute the integers $`a_n`$ via the following recursive procedure: Define $`\theta _1=\frac{1}{\theta }`$, $`a_1=\theta _1`$, $`\theta _{n+1}=\frac{1}{\theta _na_n}`$, $`a_{n+1}=\theta _{n+1}`$, where $`x`$ is the integer part of a real number $`x`$. As is customary, we use $`p_n`$, $`q_n`$ to denote the numerator and denominator, respectively, of the $`n^{th}`$ convergent:

$$[a_1,\mathrm{},a_n]=\frac{1}{a_1+\frac{1}{a_2+\frac{1}{\mathrm{}+\frac{1}{a_n}}}}=\frac{p_n}{q_n}.$$

The $`p_n`$’s and $`q_n`$’s satisfy the recursive relations

$$p_n=a_np_{n1}+p_{n2},q_n=a_nq_{n1}+q_{n2},$$

where $`p_0=0`$, $`p_1=1`$, $`q_0=1`$ and $`q_1=a_1`$. The following remarkable fact can be found in any book on continued fractions.

###### Theorem 2.1.

For every $`n`$ one has $`|\theta \frac{p_n}{q_n}|<\frac{1}{q_nq_{n+1}}<(\frac{1}{q_n})^2`$.

We now record a few trivial lemmas that will be used in later estimates.

###### Lemma 2.2.

If $`F(k)`$ is the $`k^{th}`$ Fibonacci number (where $`F(0)=1=F(1)`$, $`F(2)=2`$, etc.) then for each pair of positive integers $`n`$ and $`k`$ we have $`q_{n+k}F(k)q_n`$.

###### Proof.

Use induction and the recursion formula. ∎

###### Lemma 2.3.

$`\underset{k=0}{\overset{\mathrm{}}{}}\frac{1}{F(k)}\frac{2\sqrt{5}}{\sqrt{5}1}`$.

###### Proof.

Induction shows that $`F(k)(\frac{1+\sqrt{5}}{2})^{k1}`$ and hence the result follows by taking reciprocals and applying the formula for a geometric series. ∎

###### Lemma 2.4.

$`\underset{k=0}{\overset{\mathrm{}}{}}\frac{1}{q_{n+k}}\frac{1}{q_n}\frac{2\sqrt{5}}{\sqrt{5}1}`$.

###### Proof.

Apply the two previous lemmas. ∎

We close this section with the definition of pseudospectra.

###### Definition 2.5.

For $`\epsilon >0`$, the $`\epsilon `$-pseudospectrum of $`TB(H)`$ is

$$\sigma ^{(\epsilon )}(T)=\{\lambda :(\lambda T)^1\frac{1}{\epsilon }\},$$

where $`X^1=\mathrm{}`$ if $`XB(H)`$ is a non-invertible operator. Note that the usual spectrum $`\sigma (T)`$ is contained in $`\sigma ^{(\epsilon )}(T)`$ for every $`\epsilon >0`$ and, moreover,

$$\sigma (T)=\underset{\epsilon >0}{}\sigma ^{(\epsilon )}(T).$$

See http://web.comlab.ox.ac.uk/projects/pseudospectra for more, including software for computing pseudospectra.

## 3. The Pimsner-Voiculescu Construction

In this section we state the technical aspects of the Pimsner-Voiculescu construction that we will need. We see no reason to reproduce proofs as all the details can be found in \[5, Sections VI.4 and VI.5\].

Throughout this note, $`M_k()`$ will denote the $`k\times k`$ complex matrices. Let $`0<p<q`$ be integers, $`\omega =e^{2\pi i\frac{p}{q}}`$ and define matrices $`u_{\frac{p}{q}},v_{\frac{p}{q}}M_q()`$ as follows.

$$u_{\frac{p}{q}}=\left(\begin{array}{ccccccc}0& 1& 0& 0& \mathrm{}& 0& 0\\ 0& 0& 1& 0& \mathrm{}& 0& 0\\ 0& 0& 0& 1& \mathrm{}& 0& 0\\ 0& 0& 0& 0& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& \mathrm{}& 0& 1\\ 1& 0& 0& 0& \mathrm{}& 0& 0\end{array}\right)$$

$$v_{\frac{p}{q}}=\left(\begin{array}{ccccccc}1& 0& 0& 0& \mathrm{}& 0& 0\\ 0& \omega & 0& 0& \mathrm{}& 0& 0\\ 0& 0& \omega ^2& 0& \mathrm{}& 0& 0\\ 0& 0& 0& \omega ^3& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& \mathrm{}& \omega ^{q2}& 0\\ 0& 0& 0& 0& \mathrm{}& 0& \omega ^{q1}\end{array}\right)$$

If we are given an irrational number $`\theta (0,1)`$ and we let $`\frac{p_n}{q_n}`$ denote the $`n^{th}`$ convergent (as in the previous section) then to ease notation we will let $`\theta _n=\frac{p_n}{q_n}`$, $`u_{\theta _n}=u_{\frac{p_n}{q_n}}`$ and $`v_{\theta _n}=v_{\frac{p_n}{q_n}}`$.

###### Theorem 3.1.

(Pimsner-Voiculescu) For each natural number $`n`$ there is a $``$-monomorphism $`\phi _n:M_{q_{n1}}()M_{q_n}()M_{q_n}()M_{q_{n+1}}()`$ such that

$$\phi _n(u_{\theta _{n1}}u_{\theta _n})u_{\theta _n}u_{\theta _{n+1}}<\frac{2\pi }{q_{n1}}$$

and

$$\phi _n(v_{\theta _{n1}}v_{\theta _n})v_{\theta _n}v_{\theta _{n+1}}<\frac{\pi }{q_{n1}}+\frac{4\pi }{q_n}.$$

## 4. Computing Spectra in Irrational Rotation Algebras

For an arbitrary number $`\theta (0,1)`$ there is a universal C-algebra, denoted by $`A_\theta `$, which is generated by two unitaries $`U_\theta ,V_\theta `$ subject to the commutation relation

$$U_\theta V_\theta =e^{2\pi i\theta }V_\theta U_\theta .$$

Universality means that if $`\stackrel{~}{U}_\theta ,\stackrel{~}{V}_\theta `$ are any other unitary operators satisfying the same relation then there exists a surjective $``$-homomorphism $`A_\theta C^{}(\stackrel{~}{U}_\theta ,\stackrel{~}{V}_\theta )`$ such that $`U_\theta \stackrel{~}{U}_\theta `$ and $`V_\theta \stackrel{~}{V}_\theta `$. One crucial fact which we will need is that if $`\theta `$ is irrational then $`A_\theta `$ is a simple C-algebra (i.e. has no non-trivial, closed, two-sided ideals) and hence any $``$-homomorphism will be injective in this case (cf. \[1, Theorem 1.10\]).

In the theorem below $`d_H(,)`$ denotes the Hausdorff distance between two compact subsets of the complex plane:

$$d_H(\mathrm{\Sigma },\mathrm{\Lambda })=\mathrm{max}\{\underset{\sigma \mathrm{\Sigma }}{sup}d(\sigma ,\mathrm{\Lambda }),\underset{\lambda \mathrm{\Lambda }}{sup}d(\lambda ,\mathrm{\Sigma })\},$$

where $`d(\sigma ,\mathrm{\Lambda })=inf_{\lambda \mathrm{\Lambda }}|\sigma \lambda |`$.

###### Theorem 4.1.

Let $`\theta (0,1)`$ be irrational, $`\theta _n=\frac{p_n}{q_n}`$, $`\alpha _{\pm 1},\beta _{\pm 1}`$ be four complex numbers and $`M=\mathrm{max}\{|\alpha _{\pm 1}|,|\beta _{\pm 1}|\}`$. If we let $`H_\theta =\alpha _1U_\theta +\alpha _1U_\theta ^{}+\beta _1V_\theta +\beta _1V_\theta ^{}`$ and $`h_{\theta _n}=\alpha _1u_{\theta _n}+\alpha _1u_{\theta _n}^{}+\beta _1v_{\theta _n}+\beta _1v_{\theta _n}^{}`$ then:

1. (cf. Definition 2.5) For each $`\epsilon >0`$ we have

$$\left(\sigma ^{(\epsilon )}(h_{\theta _{n1}})\sigma ^{(\epsilon )}(h_{\theta _n})\right)\sigma ^{(\epsilon +\epsilon _n)}(H_\theta )\left(\sigma ^{(\epsilon +2\epsilon _n)}(h_{\theta _{n1}})\sigma ^{(\epsilon +2\epsilon _n)}(h_{\theta _n})\right),$$

where $`\epsilon _n=204M(\frac{1}{q_{n1}}+\frac{1}{q_n})`$.
2. (Normal Case) If $`h_{\theta _{n1}},h_{\theta _n}`$ and $`H_\theta `$ all happen to be normal operators (e.g. self-adjoints) then

$$d_H(\sigma (h_{\theta _{n1}})\sigma (h_{\theta _n}),\sigma (H_\theta ))204M(\frac{1}{q_{n1}}+\frac{1}{q_n}).$$

###### Proof.

It follows from the Pimsner-Voiculescu construction that there is an AF algebra $`𝔄_\theta `$ and a sequence of $``$-monomorphisms $`\psi _n:M_{q_{n1}}()M_{q_n}()𝔄_\theta `$ with the property that

$$\psi _n(u_{\theta _{n1}}u_{\theta _n})\psi _{n+1}(u_{\theta _n}u_{\theta _{n+1}})<\frac{2\pi }{q_{n1}}$$

and

$$\psi _n(v_{\theta _{n1}}v_{\theta _n})\psi _{n+1}(v_{\theta _n}v_{\theta _{n+1}})<\frac{\pi }{q_{n1}}+\frac{4\pi }{q_n}.$$

Due to the exponential growth of the $`q_n`$’s we have that both $`\{\psi _n(u_{\theta _{n1}}u_{\theta _n})\}`$ and $`\{\psi _n(v_{\theta _{n1}}v_{\theta _n})\}`$ are Cauchy sequences of unitaries and hence converge to some unitaries in $`𝔄_\theta `$ which satisfy the same commutation relation as $`U_\theta `$ and $`V_\theta `$. By universality there exists a (necessarily faithful – by simplicity) $``$-homomorphism $`\mathrm{\Psi }:A_\theta 𝔄_\theta `$ with the property that $`\psi _n(u_{\theta _{n1}}u_{\theta _n})\mathrm{\Psi }(U_\theta )`$ and $`\psi _n(v_{\theta _{n1}}v_{\theta _n})\mathrm{\Psi }(V_\theta )`$. More importantly, we have the following estimates:

$`\mathrm{\Psi }(U_\theta )\psi _n(u_{\theta _{n1}}u_{\theta _n})`$ $``$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\psi _{n+k}(u_{\theta _{n+k1}}u_{\theta _{n+k}})\psi _{n+k+1}(u_{\theta _{n+k}}u_{\theta _{n+k+1}})`$
$``$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2\pi }{q_{n+k1}}}`$
$`=`$ $`2\pi ({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}}+{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{q_{n+k+1}}})`$
$``$ $`2\pi ({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}})+{\displaystyle \frac{4\pi \sqrt{5}}{\sqrt{5}1}}{\displaystyle \frac{1}{q_{n+1}}}.`$
Note that we applied Lemma 2.4 to get the last inequality. A similar argument shows that

$$\mathrm{\Psi }(V_\theta )\psi _n(v_{\theta _{n1}}v_{\theta _n})\frac{\pi }{q_{n1}}+\frac{5\pi }{q_n}+\frac{10\pi \sqrt{5}}{\sqrt{5}1}\frac{1}{q_{n+1}}.$$

With these estimates in hand it is immediate that

$`\mathrm{\Psi }(H_\theta )\psi _n(h_{\theta _{n1}}h_{\theta _n})`$ $``$ $`(|\alpha _1|+|\alpha _2|)\left(2\pi ({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}})+{\displaystyle \frac{4\pi \sqrt{5}}{\sqrt{5}1}}{\displaystyle \frac{1}{q_{n+1}}}\right)`$
$`+`$ $`(|\beta _1|+|\beta _2|)\left({\displaystyle \frac{\pi }{q_{n1}}}+{\displaystyle \frac{5\pi }{q_n}}+{\displaystyle \frac{10\pi \sqrt{5}}{\sqrt{5}1}}{\displaystyle \frac{1}{q_{n+1}}}\right).`$
The estimate above, though a bit ugly, is the one which should be used in practice as the constants are smaller than those claimed in the theorem. To get the cleaner statement of the theorem we first note that the recursion formula $`q_{n+1}=a_{n+1}q_n+q_{n1}`$ implies

$$\frac{1}{q_{n+1}}\frac{1}{q_{n1}}+\frac{1}{q_n},$$

and thus

$`\mathrm{\Psi }(H_\theta )\psi _n(h_{\theta _{n1}}h_{\theta _n})`$ $``$ $`2M\left(2\pi +{\displaystyle \frac{4\pi \sqrt{5}}{\sqrt{5}1}}\right)({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}})`$
$`+`$ $`2M\left(5\pi +{\displaystyle \frac{10\pi \sqrt{5}}{\sqrt{5}1}}\right)({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}})`$
$`=`$ $`2M\left(7\pi +{\displaystyle \frac{14\pi \sqrt{5}}{\sqrt{5}1}}\right)({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}})`$
$`=`$ $`14\pi M{\displaystyle \frac{3\sqrt{5}1}{\sqrt{5}1}}({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}})`$
$``$ $`204M({\displaystyle \frac{1}{q_{n1}}}+{\displaystyle \frac{1}{q_n}}).`$
Now we quote some general results concerning approximation of spectral quantities. Given two operators $`S,TB(H)`$ one has:

1. (cf. \[6, Theorem 3.27\]) For every $`\epsilon >0`$, $`\sigma ^{(\epsilon )}(S)\sigma ^{(\epsilon +ST)}(T)\sigma ^{(\epsilon +2ST)}(S)`$.
2. If both $`S`$ and $`T`$ happen to be normal then $`d_H(\sigma (S),\sigma (T))ST`$.

With these general facts and the norm estimates above one easily completes the proof by observing that the spectral quantities of $`ST`$ (i.e. pseudo or actual spectrum) are just the union of the spectral quantities for $`S`$ and $`T`$. ∎

###### Remark 4.2.

Almost tri-diagonality. Note that $`h_{\theta _n}`$ is “almost” a tri-diagonal matrix. It is non-zero in the upper-right and lower-left entries but otherwise is tri-diagonal. In particular, these matrix models are sparse.

###### Remark 4.3.

Discretized Schrödinger operators with polynomial potential.

The theorem above is stated for a class of operators that includes all of the examples discussed in  (i.e. Almost Mathieu operators and a few non-self-adjoint examples). However, similar results hold for a broader class of examples that includes all the one-dimensional discretized Schrödinger operators with polynomial potential. In fact, if $`R(X,Y)`$ is any polynomial in non-commuting variables $`X`$ and $`Y`$ then

$$\mathrm{\Psi }(R(U_\theta ,V_\theta ))\psi _n(R(u_{\theta _{n1}},v_{\theta _{n1}})R(u_{\theta _n},v_{\theta _n}))=O(\frac{1}{q_{n1}}+\frac{1}{q_n}).$$

Unfortunately, the constants involved get quite complicated quite quickly but they can be explicitly worked out for a given polynomial should one so desire. Having this norm estimate one gets similar convergence results for spectral quantities in this setting. Finally we should point out that the spectral quantities of any operator in an irrational rotation algebra (e.g. Schrödinger operators with arbitrary quasiperiodic potential) can, in principle, be approximated with this method – first approximate by a polynomial and then proceed – however, it seems impossible to control rates of convergence in this generality.

###### Remark 4.4.

These rates of convergence are sharp.

To see this we consider the case that $`H_\theta =U_\theta `$ and $`\theta `$ is an irrational number for which the $`a_n`$’s are uniformly bounded. Then $`d_H(\sigma (H_\theta ),\sigma (h_{\theta _{n1}})\sigma (h_{\theta _n}))292(\frac{1}{q_{n1}}+\frac{1}{q_n})`$. It is easy to check that $`\frac{1}{q_{n1}}+\frac{1}{q_n}\frac{4+a_n}{q_{n1}+q_n}`$ and hence if one could prove that $`d_H(\sigma (H_\theta ),\sigma (h_{\theta _{n1}})\sigma (h_{\theta _n}))=o(\frac{1}{q_{n1}}+\frac{1}{q_n})`$ then it would follow that $`d_H(\sigma (H_\theta ),\sigma (h_{\theta _{n1}})\sigma (h_{\theta _n}))=o(\frac{1}{q_{n1}+q_n})`$. However this is impossible as it would imply that the Lebesgue measure of $`\sigma (H_\theta )`$ is zero (which it isn’t: $`\sigma (U_\theta )=\{\lambda :|\lambda |=1\}`$) since there are at most $`q_{n1}+q_n`$ points in $`\sigma (H_{\theta _{n1}})\sigma (H_{\theta _n})`$.

###### Remark 4.5.

Implementability. It is, of course, true that for some irrational numbers the strategy proposed in this paper for approximating spectra can’t reasonably be implemented. For example, it can happen that $`q_1`$ is very small while $`q_2`$ is far larger than any computer can handle (just take $`a_2`$ to be huge). On the other hand, the set of irrationals for which the associated sequences $`\{a_n\}`$ stay bounded is a dense set in $`(0,1)`$. For example, all irrationals which are the roots of quadratic equations have this property; hence $`\sqrt{r}`$ for any rational number and these are already dense in $`(0,1)`$. In other words, for a dense set of irrationals this program is reasonable and hence the question becomes how spectra behave as the parameter $`\theta `$ changes. But it is well known that spectra of operators in irrational rotation algebras vary continuously in $`\theta `$ (see the next section) and hence the program suggested in this paper is about as practical as could be hoped for.

###### Remark 4.6.

Uniqueness of convergents. Another basic, but remarkable, property of continued fractions is that if $`\theta `$ is irrational and $`|\theta \frac{p}{q}|<\frac{1}{2q^2}`$ then $`\frac{p}{q}`$ is necessarily one of the convergents in the continued fraction decomposition of $`\theta `$ (i.e. if $`\theta =[a_1,a_2,\mathrm{}]`$ then there is an $`n`$ such that $`\frac{p}{q}=[a_1,\mathrm{},a_n]=\frac{p_n}{q_n}`$ – see \[7, Theorem 184\]). The point is that if one starts with a rational $`\frac{p}{q}`$ and only looks at the irrationals $`\theta `$ which are a distance at most $`\frac{1}{2q^2}`$ away from $`\frac{p}{q}`$ then the spectral quantities of $`h_{\frac{p}{q}}`$ will provide good approximations to those of $`H_\theta `$ (since $`\frac{p}{q}`$ necessarily arises in the continued fraction decomposition of all such $`\theta `$ and hence the theorem above applies). Of course, rates of convergence become trickier in this setting but our only point is that one can simultaneously approximate a small interval of operators with a single matrix model.

## 5. General One-sided Rates of Convergence

In our final section we will record a consequence of the following theorem of Haagerup and Rørdam (cf. \[1, Corollary 3.5\]). Among other things this shows that if one fixes a polynomial in the canonical generators then the spectral quantities vary continuously in $`\theta `$.

###### Theorem 5.1.

For every $`\theta ,\theta ^{}[0,1)`$ there exists a Hilbert space $`H`$ and injective $``$-homomorphisms $`\pi _\theta :A_\theta B(H)`$, $`\pi _\theta ^{}:A_\theta ^{}B(H)`$ such that

$$\pi _\theta ^{}(U_\theta ^{})\pi _\theta (U_\theta ),\pi _\theta ^{}(V_\theta ^{})\pi _\theta (V_\theta )9\sqrt{6\pi |\theta \theta ^{}|}.$$

As in the last section we will restrict attention to operators of the form $`H_\theta =\alpha _1U_\theta +\alpha _1U_\theta ^{}+\beta _1V_\theta +\beta _1V_\theta ^{}`$ as the estimates are cleaner in this setting. However, the same techniques handle more general polynomials as well. In the theorem below we only get “one-sided” rates of convergence (as opposed to estimates on the Hausdorff distance). The modest advantage, however, is that this strategy is implementable for any irrational number.

We find the following notation convenient: For compacts $`\mathrm{\Lambda },\mathrm{\Sigma }`$ and $`\delta >0`$ we write $`\mathrm{\Lambda }^\delta \mathrm{\Sigma }`$ if for each $`\lambda \mathrm{\Lambda }`$ there exists $`\sigma \mathrm{\Sigma }`$ such that $`|\lambda \sigma |<\delta `$. Note that $`d_H(\mathrm{\Lambda },\mathrm{\Sigma })<\delta `$ if and only if $`\mathrm{\Lambda }^\delta \mathrm{\Sigma }`$ and $`\mathrm{\Sigma }^\delta \mathrm{\Lambda }`$.

###### Theorem 5.2.

Let $`\theta (0,1)`$ be irrational and for each $`n`$ choose an integer $`p_n\{0,1,\mathrm{},n1\}`$ such that $`|\theta \frac{p_n}{n}|\frac{1}{2n}`$. If $`H_\theta `$ is an operator polynomial as above, $`h_n=\alpha _1u_{\frac{p_n}{n}}+\alpha _1u_{\frac{p_n}{n}}^{}+\beta _1v_{\frac{p_n}{n}}+\beta _1v_{\frac{p_n}{n}}^{}`$ is the associated matrix model and $`M=\mathrm{max}\{|\alpha _{\pm 1}|,|\beta _{\pm 1}|\}`$ then the following statements hold:

1. For every $`\epsilon >0`$,

$$\sigma ^{(\epsilon )}(h_n)\sigma ^{(\epsilon +\frac{C_1}{\sqrt{n}})}(H_\theta ),$$

where $`C_136M\sqrt{3\pi }.`$
2. If $`H_\theta `$ happens to be normal then

$$\sigma (h_n)^{\frac{C_1}{\sqrt{n}}}\sigma (H_\theta ),$$

where $`C_136M\sqrt{3\pi }.`$

###### Proof.

To prove this result we first observe that all of the spectral quantities associated to $`h_n`$ are contained in the corresponding quantities for $`H_{\frac{p_n}{n}}=\alpha _1U_{\frac{p_n}{n}}+\alpha _1U_{\frac{p_n}{n}}^{}+\beta _1V_{\frac{p_n}{n}}+\beta _1V_{\frac{p_n}{n}}^{}`$ (since universality gives us a $``$-homomorphism taking $`H_{\frac{p_n}{n}}h_n`$). The Lip<sup>1/2</sup>-continuity theorem of Haagerup-Rørdam then gives upper bounds on the distance between the spectral quantities of $`H_\theta `$ and $`H_{\frac{p_n}{n}}`$ which completes the proof. ∎

###### Remark 5.3.

Hausdorff Convergence. It can be shown that $`d_H(\sigma ^{(\epsilon )}(h_n),\sigma ^{(\epsilon )}(H_\theta ))0`$ (for every $`\epsilon >0`$). The proof is similar to the proof of \[4, Theorem 3.5\] and hence will be omitted. However, it does not seem possible to control the rate of convergence in the Hausdorff metric. For example, it is certainly not the case that rates of convergence are controlled by $`|\theta \frac{p_n}{n}|^\sigma `$ for any $`\sigma >0`$ because many irrational numbers (e.g. transcendental numbers) have extremely good rational approximations. To be more precise, it can happen that a particular $`\theta `$ has the property that there are infinitely many co-prime solutions to the inequality $`|\theta \frac{p}{q}|<\frac{1}{q^k}`$ where $`k`$ is some fixed natural number. If the Hausdorff distance between spectra could be bounded above by $`|\theta \frac{p}{q}|^\sigma `$ then we would also get an upper bound of $`(\frac{1}{q})^{k\sigma }`$. However, when $`k\sigma >1`$ it would then follow that the spectrum of $`H_\theta `$ has Lebesgue measure zero since there are at most $`q`$ eigenvalues for the matrix $`h_{\frac{p}{q}}`$. Since these operators need not have spectra of measure zero it follows that an upper bound of the form $`|\theta \frac{p_n}{n}|^\sigma `$ is impossible.


format: text

text:
# Evidence for a 𝑘^{-5/3} spectrum from the EOLE Lagrangian balloons in the low stratosphere

## I Introduction

The EOLE project consisted in the release of 483 constant-volume pressurized balloons, in the Southern Hemisphere mid-latitudes throughout the period September 1971-March 1972, at approximately 200 hPa. The life-time of these balloons was from a few days to about one year, with an average value of about 120 days. Their motion was basically isopycnal except for small diurnal volume variations of the envelop of less than 1% and inertial oscillations of a few meters in the vertical, excited by wind shear and small-scale turbulence. The position of the balloons and meteorological parameters were periodically transmitted to satellite by ARGOS system. The trajectories of the EOLE experiment still provide nowadays the most extensive data set of experimental quasi-Lagrangian tracers in the atmosphere for observing the properties of the medium-to-large-scale motion at the tropopause level.

Both Eulerian and Lagrangian analyses have been performed by several authors. Morel and Desbois (1974) deduced the mean circulation around 200 hPa from the balloon flights, as formed by a mid-latitude zonal stream with a meridional profile characterized by a typical velocity $`30`$ $`\mathrm{ms}^1`$ inside the jet, overimposed to meridional velocity field disturbances of much smaller intensity, $`1`$ $`\mathrm{ms}^1`$, and to residual standing waves acting as spatial perturbations of the zonal velocity pattern, producing the typical shape of a meandering jet. These results have been largely confirmed by operational analysis since then. Morel and Larchevêque (1974), hereafter ML, investigated the synoptic-scale turbulent properties. They measured the mean square relative velocity and the relative diffusivity of balloon pairs, and found essentially two major regimes for Lagrangian dispersion: exponential separation for time delays less than 6 days and standard diffusion for larger times. These authors also observed that the scaling of the relative diffusivity with the separation length between two balloons agreed with a direct 2D turbulent cascade, with energy spectrum $`E(k)k^3`$, or steeper, in the range 100-1000 km. Further Eulerian analyses of large-scale velocity spectra by Desbois (1975) were compatible with the scenario proposed by Morel and Larchevêque (1974) about isotropic and homogeneous 2D turbulence with a $`k^3`$ energy distribution up to scales $`1000`$ km.

Later, other authors reached for different conclusions after observing energy spectra in the low stratosphere, measured from experimental data recorded from commercial aircraft flights, Gage (1979), Lilly (1983), Nastrom and Gage (1985). Their picture suggested a 2D turbulent inverse cascade, characterized by the $`E(k)k^{5/3}`$ spectrum, inside the interval of scales 100-1000 km.

More recently, Lindborg and Cho (2000, 2001a and 2001b) computed velocity spectra using data recorded during the MOZAIC program and also found a $`k^{5/3}`$ spectrum. They suggested a dynamical mechanism different from 2D inverse cascade where energy is injected in the large scales by breaking Rossby-gravity waves and generates a chain process down to smaller scales. Their hypothesis is supported by the observation of a downscale energy flux, whereas 2D inverse cascade should exhibit upscale energy flux (Lindborg and Cho, 2000).

We wanted to reconsider this issue by performing a new analysis of the relative dispersion properties of the EOLE balloons within the framework of dynamical system theory. Relative dispersion properties are analyzed through the computation of the Finite-Scale Lyapunov Exponent, or shortly FSLE (Aurell et al. 1997, Artale et al. 1997, Boffetta et al. 2000a). The FSLE is based on the growth rate statistics of the distance between trajectories at fixed scale, and is a better tool at analyzing scale dependent properties than plain dispersion, as explained below. This new method has been already exploited for studies of relative dispersion in atmospheric and oceanic systems (Lacorata et al. 2001, Joseph and Legras 2002, Boffetta et al. 2001, LaCasce and Ohlmann 2003, Gioia et al. 2003) and also in laboratory convection experiments (Boffetta et al. 2000b).

This paper is organized as follows: in section II we describe the FSLE methodology; section III contains the results obtained from our analysis of the EOLE experimental data; in section IV we discuss the physical information that can be extracted from this paper and possible perspectives.

## II Finite-scale relative dispersion

Generally speaking, most flows exhibit a range of scales over which fluid motion is expected to display different characteristics: a small-scale range where the velocity variations can be considered as a smooth function of space variables; a range of intermediate lengths corresponding to the coherent structures (and/or spatial fluctuations) present in the velocity field over which velocity variations are rough but highly correlated; a large-scale range over which spatial correlations have definitely decayed. In each of these ranges, relative dispersion between trajectories is governed by a different physical mechanism (chaos, turbulence, diffusion) which can be, in principle, identified from the observations. In fully developed three-dimensional turbulence, motion is only smooth under the Kolmogorov dissipative scale. In the free stratified atmosphere (above the planetary boundary layer), turbulence is a relatively rare event: motion is most often smooth but for some localized convective or turbulent events, associated with mesoscale systems, that mix momentum and tracers. Hence one expects to find a smooth (chaotic) dispersion range ended at a scale characteristic of the spacing of mixing events, followed by a range covering the mesoscale to synoptic range, and finally standard diffusion at planetary scale. This view is supported by the ubiquitous observation of long-lived laminated structures in the free troposphere (Newell et al., 1999).

In order to fix some terminology, we will use both symbols $`R`$ and $`\delta `$ for indicating the distance between balloons: the former will be considered as a quantity function of time, the latter as an argument for scale-dependent functions.

### II.1 Diffusion and Chaos

Diffusion is characterized in terms of diffusion coefficients, related to the elements of a diffusion tensor defined as

$$D_{ij}=\underset{t\mathrm{}}{lim}\frac{1}{2t}(x_i(t)x_i(0))(x_j(t)x_j(0))$$
(II.1)
where $`x_i(t)`$ and $`x_j(t)`$ are the $`i`$-th and $`j`$-th coordinates at time $`t`$, with $`i,j=1,2,3`$. The average operation $``$ is meant to be performed over a large number of particles. The diagonal elements are just the diffusion coefficients. When the $`D_{ii}`$’s are finite, then the diffusion is standard. This means that at long times, after the Lagrangian velocity correlations have decayed (Taylor, 1921), the variance of the particle displacement follows the law:

$$𝐱(t)𝐱(0)^22Dt$$
(II.2)
In presence of a velocity field characterized by coherent structures, it is more useful to observe the relative dispersion between the trajectories, rather than the absolute dispersion from the initial positions, given by (II.2), which is unavoidably dominated by the mean advection.

In the case of the EOLE experiment, where observing the expansion of balloon clusters with more than two elements is a rare event (ML) a measure of relative dispersion is given by the mean square inter-particle distance:

$$R^2(t)=𝐱^{(m)}(t)𝐱^{(n)}(t)^2$$
(II.3)
averaged over all the pairs $`(𝐱^{(m)},𝐱^{(n)})`$, where $`m`$ and $`n`$ label all the available $`N`$ trajectories. Notice that the norm in (II.2) and (II.3) must be defined accordingly to the geometry of the fluid domain, i.e. in the atmosphere we use the arc distance on the great circle of the Earth connecting the two points. The quantity in (II.3) can be measured for both initially close pairs, balloons released from the same place at short time delay, and so-called chance pairs, balloons initially distant which come close to each other at a certain time and, later, spread away again (ML). Consistency of the average in eq. (II.3) requires all the trajectory pairs to have nearly the same initial distance, a condition which strongly limits the statistics. At long times, $`R(t)^2`$ defined in eq. (II.3) is expected to approach the function $`4Dt`$, where the $`4`$ factor accounts for relative diffusion. When it happens that $`R(t)^2t^{2\nu }`$ with $`\nu >1/2`$, instead, the Lagrangian dispersion is considered as super-diffusion. A well-known example is the Richardson’s law for the particle pair separation in 3D turbulence, for which $`\nu =3/2`$ (Richardson 1926; Monin and Yanglom 1975).

On the other hand, in the limit of infinitesimal trajectory perturbations, much smaller than the characteristic lengths of the system, the evolution of the particle pair separation is characterized by the Lyapunov exponent (Lichtenberg and Lieberman, 1982), such that

$$\lambda =\underset{t\mathrm{}}{lim}\underset{R(0)0}{lim}\frac{1}{t}\mathrm{ln}\frac{R(t)}{R(0)}$$
(II.4)
If $`\lambda >0`$ the growth is exponential and the motion is said chaotic. Chaos is a familiar manifestation of non linear dynamics, leading to strong stirring of trajectories (Ottino, 1989). The process, for example, of repeated filamention around the polar vortex is basically due to Hamiltonian chaos (Legras and Dritschel, 1993). For finite perturbations within a smooth flow, the properties of exponential separation are observed for a finite time.

### II.2 Finite-Scale Lyapunov Exponent

The idea of FSLE (Aurell et al. 1997; Artale et al. 1997), was formerly introduced in the framework of the dynamical systems theory, in order to characterize the growth of non-infinitesimal perturbations (i.e. the distance between trajectories). If $`\delta `$ is the scale of the perturbation, and $`\tau (\delta )`$ is the mean time that $`\delta `$ takes to grow a factor $`r>1`$, then the FSLE is defined as

$$\lambda (\delta )=\frac{1}{\tau (\delta )}\mathrm{ln}r$$
(II.5)
The average operation is assumed to be performed over a large ensemble of realizations. For factors $`r`$ not much larger than $`1`$, $`\lambda (\delta )`$ does not depend sensitively on $`r`$. If $`r=2`$ then $`\tau (\delta )`$ is also called doubling time. Operatively, $`N+1`$ scales are chosen to sample the spatial range of perturbation, $`\delta _0<\delta _1<\mathrm{}<\delta _N`$, and a growth factor $`r`$ is defined such that $`\delta _i=r\delta _{i1}`$ for $`i=1,N`$. Let $`l_{min}`$ and $`l_{max}`$ be the smallest and the largest characteristic length of the system, respectively. If $`\delta _0l_{min}`$ then the FSLE characterizes the doubling time of infinitesimal perturbation. In the opposite side of the range, if $`\delta _Nl_{max}`$ then the FSLE follows the scaling law of diffusion $`\lambda (\delta )\delta ^2`$ for $`\delta \delta _N`$, as can be deduced by noticing that the mean square particle distance must grow linearly in time, see (II.2). In general, if the mean square size of a tracer concentration follows the $`R^2t^{2\nu }`$ law, the FSLE scales as $`\lambda (\delta )\delta ^{1/\nu }`$. As we have seen before, for standard diffusion $`\nu =1/2`$ while for Richardson’s super-diffusion $`\nu =3/2`$. The main interest of FSLE is to select processes occurring at a fixed scale. We stress that definition (II.5) differs substantially from

$$\lambda ^{}(\delta )=\frac{1}{R^2}\frac{dR^2}{dt}|_{R^2=\delta ^2}$$
(II.6)
defined in terms of the mean square relative displacement, because of the different averaging procedures in the two cases: $`R^2`$ is computed at fixed time while $`\tau (\delta )`$ is computed at fixed scale. As a result, a physical situation which is well characterized in terms of FSLE, either for scaling properties or the existence of transport barriers, may be less easily characterized by studying the time growth of trajectory separation (Boffetta et al., 2000b; Joseph and Legras, 2002). One reason is that $`R^2(t)`$ depends on contribution from different regimes, as seen, for example in $`3D`$ turbulence where a dramatic dependence of $`R^2(t)`$ upon $`R^2(0)`$ is observed, even at very large Reynolds number (Fung et al., 1992)

In cases where advection is strongly anisotropic, e.g. in presence of a structure like the stratospheric jet stream, it may be useful to define the FSLE in terms of meridional (cross-stream) displacement only:

$$\lambda _{mer}(\delta ^{(mer)})=\frac{1}{\tau (\delta ^{(mer)})}\mathrm{ln}r$$
(II.7)
where $`\delta ^{(mer)}`$ is the latitude distance (or meridian arc) between two points.

Informations about the relative dispersion properties are also extracted by another fixed-scale statistics, the Finite-Scale Relative Velocity (FSRV), named by analogy with FSLE, that is defined as

$$\nu _2(\delta )=\delta 𝐯(\delta )^2$$
(II.8)
where

$$\delta 𝐯(\delta )^2=(\dot{𝐱}^{(1)}\dot{𝐱}^{(2)})^2$$
(II.9)
is the square Lagrangian velocity difference between two trajectories, $`𝐱^{(1)}`$ and $`𝐱^{(2)}`$, on scale $`\delta `$, that is for $`|𝐱^{(1)}𝐱^{(2)}|=\delta `$. The FSRV can be regarded as the 2nd order structure function of the Lagrangian velocity difference and provides a complementary analysis to the FSLE diagnostics. In particular, in the regime of Richardson’s super-diffusion, the expected behavior for the FSRV is $`\nu _2(\delta )\delta ^{2/3}`$.

We report in the next section the results of our analysis.

## III Analysis of the EOLE Lagrangian data

After a preliminary data check, the number of balloons selected for the analysis has been reduced to 382. This has been obtained by discarding ambiguous ident numbers (some ident numbers have been used twice during the campaign), discarding trajectories that cross the equator and short tracks of less than 10 points.

Successive points along a balloon trajectory were mostly recorded at a time interval of $`10^1`$ day ($`2.4`$ hours), but the overall distribution of the raw data does not cover uniformly the time axis. Hence, each of the coordinates (longitude and latitude) of every balloon trajectory has been interpolated in time by a cubic spline scheme, with a sampling rate of 25 points per day. Because of possible data impurities, each Lagrangian velocity value is monitored at every time step ($`0.04`$ day) and data segment with abnormally fast motions are discarded.

As pointed out by ML, a way to measure the dispersion between balloons is waiting for one of them to get close to another one, at a distance less than a threshold $`\delta _0`$, and then observing the evolution of their relative distance in time. This procedure is repeated for each balloon trajectory until the whole set of pairs is analyzed. The dataset includes original pairs of balloons that were launched within a short time interval and chance pairs of balloons meeting suddenly after a number of days. For the largest values of the threshold used in this study, the number of chance pairs largely exceeds the number of original pairs. In this way, global properties of the Lagrangian transport are extracted from the contributions of balloon pairs randomly distributed all over the domain. The number of balloon pairs and its evolution as the separation crosses the $`N`$ scales defined above is described by Table 1.

In Figure 1 four global relative dispersion curves are plotted, referring to four different initial thresholds $`\delta _0=`$ 25 km, 50 km, 100 km and 200 km. The statistical samples vary roughly in proportional way with $`\delta _0`$. Relative dispersion depends sensitively, as expected, on the initial conditions; the four curves meet together for separation larger than about 2500 km and saturation begins for separation larger than 4500 km, leaving room only for a short standard diffusive regime between these two separations and over a time duration of less than 10 days. The eddy diffusion coefficient, $`D_E`$, estimated by fitting the linear law $`4D_Et`$, results in $`D_E\mathrm{2.9\hspace{0.33em}10}^6`$ m2s-1, a value compatible with what was found by ML. The pre-diffusive regime is not very clear, we can say that the behavior of the balloon separation looks like a power law with exponent (changing in time) between 3 and 1.

We report in Figure 2 the mean logarithmic growth of the balloon relative separation over all pairs selected by the threshold 25 km. At very short times ($`<`$ 1 day) the slope corresponds to an exponential growth rate with $`e`$-folding time $``$ 0.4 day that we consider as a rough estimate of the LLE. At later times, the slope gradually decreases as the separation growth tends to a power law regime. In the same figure we also show the mean logarithmic growth of the inter-balloon distance computed for two 4-element clusters (that we label as ’A’ and ’B’), launched with a time interval of 3 days between them. A linear behavior (exponential growth) for both clusters is observed for short intervals; we observe as the exponential regime lasts longer for the ’A’ cluster ($`3`$ days) than for the ’B’ cluster ($`1`$ days). This illustrates the fact that the duration of a dispersion regime, here the chaotic one, may exhibit large fluctuations generally due to different meteorological conditions. As a result, average time-dependent quantities, like $`R^2(t)`$, sample different regimes at once and are poor diagnostics of dispersion properties. Incidentally, in ML the behavior of the relative dispersion between 100 and 1000 km is fitted by means of an exponential with characteristic e-folding time $`2.7`$ days (see ML, figure 8), which is compatible with the growth rate of Figure 2 between the two horizontal lines marking the range 100-1000 km, if one wants to fit it with an exponential curve branch.

Figure 3 shows the global FSLE relatively to the same four initial thresholds used for the relative dispersion and setting the amplification ratio $`r`$ to $`\sqrt{2}`$. The main result of this study is that up to about 1000 km there is evidence of Richardson’s super-diffusion, compatible with a $`k^{5/3}`$ spectrum, displayed by the behavior $`\lambda (\delta )=\alpha \delta ^{2/3}`$. The best fit is obtained for the initial thresholds 100 and 200 km which encompass a much larger number of pairs than smaller thresholds (see Table 1). The quantity $`\alpha ^3`$ is physically related to the mean relative kinetic energy growth rate (for unit mass) between balloons moving apart. Standard diffusion is approached at scales larger than 2000 km. The value of the eddy diffusion coefficient is estimated by fitting the FSLE in the diffusive range with $`(4\mathrm{ln}r)D_E\delta ^2`$, as shown in Boffetta et al. (2000a) by means of a dimensional argument. We find that this value is $`D_E10^7`$ m2s-1. Notice that the initial threshold does not affect very much the general behavior, except for obvious changes in the statistical samples.

Figure 4 shows global (mainly zonal) and meridional ($`\lambda _{mer}`$, see (II.7)) FSLE of the balloon pairs with initial threshold $`100`$ km. We find that the dispersion is basically isotropic up to scales of about 500 km, which is in rough agreement with the results of Morel and Larchevêque (they give a value three times larger but their analysis, see their Fig. 7, does not display a well-defined cut-off). At scales larger than 500 km, the two components of the FSLE decouple and the meridional dispersion rate follows the standard diffusion law $`\delta ^2`$ with a meridional eddy diffusion coefficient $`D_E10^6`$ m2s-1.

In order to compute the FSRV, the relative velocity between balloons is approximated by the finite difference formula $`(|R(t+\mathrm{\Delta }t)||R(t)|)/\mathrm{\Delta }t`$, where $`|R(t)|`$ is the absolute value of the great circle arc between two balloons at time $`t`$ and $`\mathrm{\Delta }t=0.04`$ day is the time interval between two successive points along a trajectory. The properties of the Lagrangian relative velocity are shown in Figure 5. The FSRV confirms the results obtained with the FSLE: between 100 and 1000 km the behavior is $`\delta ^{2/3}`$, corresponding to the Richardson’s law; asymptotic saturation sets in beyond this range (fully uncorrelated velocities).

## IV Discussion and conclusions

We have revisited the dispersion properties of the EOLE Lagrangian dataset using a new approach, using Finite-Scale Lyapunov Exponent, that is better suited to analyze scale dependent properties than standard tools that were used, e.g., by ML in a previous study of the same dataset. We were motivated by the fact that ML found results supporting a $`k^3`$ inertial range between 100 and 1000 km, whereas more recent studies based on aircraft data found a $`k^{5/3}`$ behavior in the same range of scales.

Our main result of our improved analysis is that the EOLE dataset supports a $`k^{5/3}`$ behavior in the range 100-1000 km as shown by the scaling properties of FSLE in this range indicating Richardson’s superdiffusion. At distances smaller than 100 km, our results suggest an exponential separation with an e-folding time of about one day, in rough agreement with ML. At scales larger 1000 km, the dispersion tends to a standard diffusion before saturating at the planetary scale. Since the large-scale flow is dominated by the meandering zonal circulation, estimated diffusion coefficient is 10 times larger for total dispersion ($`10^7`$ m2s-1) than for meridional dispersion ($`10^6`$ m2s-1).

Our result is compatible with an inverse 2D energy cascade in the range 100-1000 km or with the recently proposed alternative of a direct energy cascade (Lindborg and Cho, 2000). Our study of the EOLE experiment has shown that this still unparalleled dataset of Lagrangian trajectories in the atmosphere is in agreement with results obtained from aircraft data. The challenge is now to compare these trajectories with the global wind fields produced by the recent reanalysis by operational weather centers.

###### Acknowledgements.

We warmly thank G. Boffetta, F. D’Andrea, V. Daniel, B. Joseph, A. Mazzino, F. Vial for interesting discussions and suggestions. G.L. thanks the European Science Foundation for financial support through a TAO Exchange Grant 1999, and the LMD-ENS (Paris) for hosting him.

FIGURE CAPTIONS

Figure 1. Mean square balloon separation. The four curves refers to 4 different initial thresholds: 25 km (a), 50 km (b), 100 km (c) and 200 km (d). All the curves but for 25 km have been shifted in time in order to collapse together for $`R^2>10^7`$ km2. The eddy diffusion coefficient corresponding to the indicated slope is $`D_E\mathrm{2.9\hspace{0.33em}10}^6`$ m2s-1. Units: time in days, $`R^2`$ in km2.

Figure 2. Mean logarithmic growth of the balloon separation ($`+`$). The initial separation is $``$ 25 km. The two clusters ’A’ ($``$) and ’B’ ($`\times `$) have 4 balloons each and were launched with a 3-day lag in November 1971. The $`e`$-folding time of the exponential growth is $`0.4`$ day. The two horizontal lines mark the range 100-1000 km (natural logarithm units). The straight line ML is the result found by Morel and Larchvêque (1974) in their fig. 8. Units: time in days, $`R`$ in km.

Figure 3. FSLE of the balloon pairs, four curves with the same four initial thresholds as in Figure 1. The eddy diffusion coefficient is $`D_E10^7`$ m2s-1. The quantity $`\alpha ^3`$ gives the order of magnitude of the relative kinetic energy growth rate (for unit mass) between balloons in the Richardson’s regime $`\delta ^{2/3}`$. Units: $`\delta `$ in km, $`\lambda `$ in day-1.

Figure 4. FSLE of the balloon pairs, describing total ($``$) and meridional ($`\times `$) dispersion, with initial threshold $`100`$ km. The meridional FSLE is $`\lambda _{mer}`$ defined in (II.7). The meridional eddy diffusion coefficient is $`D_E\mathrm{1.5\hspace{0.17em}10}^6`$ m2s-1. Units: $`\delta `$ in km, $`\lambda `$ in day-1.

Figure 5. FSRV of the balloon pairs for initial threshold 50 km. The reference velocity is $`u_0=`$100 km/hour. The slope $`2/3`$ corresponds to the Richardson’s law. Units: $`\delta `$ in km.


text:
# Specific protein-protein binding in many-component mixtures of proteins

## 1 Introduction

Proteins must interact to function, and they must interact with very high specificity in the highly complex mixtures that lie inside cells. For example, the 3 proteins Cdk2, cyclin and PCNA bind together in the nucleus to form a complex, which then goes on to interact with high specificity with a few other proteins . (PCNA=Proliferating Cell Nuclear Antigen and Cdk2=Cyclin dependent kinase; they are both involved in cell division but also in other processes.) Cdk2 must bind strongly to cyclin and to PCNA, i.e., the complex must have a large equilibrium constant or equivalently a small dissociation constant. But Cdk2 must not bind to the of order one thousand other proteins present i.e., the equilibrium constants with all these other proteins must be small. This highly specific and strong binding in a mixture of thousands of different molecules is a demanding requirement and here we will try to quantify how demanding it is for a very simple model of a protein. To do so, we will rely on the fact that specific binding in a mixture is analogous to communication through a noisy channel, which has been extensively studied, in particular by Shannon  and the many who have built on his work. In communication through a noisy channel reliable communication requires that all possible messages be sufficiently distinct from each other even after distortion by the noise. In specific binding in a mixture, a protein binding site or patch on the surface of a protein must be sufficiently distinct from all other patches to allow one other patch to bind to it with an equilibrium constant much larger than that of any other interaction of the patch. The theory developed by Shannon has been enormously useful in understanding communication systems and building efficient ones. We hope that applying the theory to protein-protein interactions in cells will be similarly useful.

This analogy to communication has been used for the similar problem of binding between a protein and a specific DNA sequence, in particular by Schneider . See also  for other related work. Similar approaches have also been used extensively, see the review , in studies of our immune system, see  for an introduction to our immune system. Studies of the immune system look at a large set of immune-system proteins recognising a single foreign protein. By contrast, our situation, a model of the cytoplasm of a prokaryote or a compartment of a eukaryote, involves not a single protein binding to another but binding between a very large number of pairs of proteins. However, there are close analogies between the specific binding required inside all cells and that required of immune-system proteins, in particular in both cases proteins must avoid binding to the wrong protein . Inside cells, this would be all proteins bar a protein’s partner, whereas in the case of the immune system, its proteins must not bind to other human proteins.

Our example of Cdk2 needing to bind to PCNA and to cyclin, is the rule rather than the exception in eukaryote cells . We merely chose a specific example to make the discussion more concrete. The cells rely totally on a complex network or web of many specific interactions. This large set of interactions is often called the interactome, by analogy to the use of the word genome to denote the set of genes of an organism. In the study of noisy channels, a fundamental quantity of interest is the channel capacity : the maximum number of bits of information that can be sent per second. Shannon derived his famous channel-capacity theorem  which allows this channel capacity to be calculated. Here we will calculate a related quantity for proteins, which we will call the capability. It is the maximum number of specific interactions, the maximum size of the interactome.

Protein interactions need to be understood to understand how cells work, but also the specificity of interactions is crucial to drug development. For example, Colas et al.  used a combinatorial-library approach to obtain peptide aptamers that bind specifically to Cdk2. Peptide aptamers are proteins that consist of a basic globular protein to form a ‘scaffold’ plus a region where the amino acids are allowed to vary. A library of these is generated and the peptide aptamers whose variable regions bind most strongly to the selected protein, here Cdk2, are selected. To function effectively, the potential drug, here the peptide aptamer, must interact in a highly specific way. If the requirement is to interfere with Cdk2, then the drug should bind to Cdk2 but not to the other thousands of proteins present. It is not sufficient that it bind strongly to Cdk2 in a dilute solution in vitro.

In the next section we will specify our mixture that is a model for the cytoplasm. It is an $`N`$ component mixture with highly specific interactions between the components. Section 3 defines the model of a protein. It is a simple lattice model in which each model protein has 6 binding sites or surface patches that mediate the specific interactions. Each patch is specified by a sequence of elements, and the interactions of a pair of patches is determined by their sequences. The model is not new, it was introduced in , but there only solubility, not specific binding was considered. The capability, the maximum number of interactions of a given high specificity, is calculated in section 4, as a function of the length of these sequences. The next section, section 5, contains some illustrative results. The final section is a conclusion. An appendix contains a discussion of Shannon’s channel-capacity theorem and a brief analysis along the lines of Schneider .

## 2 Model of the cytoplasm of a cell

Our picture of the mixture of proteins inside a cell is of a large number of proteins, interacting via interactions of high specificity . Our model of this is a mixture with $`N`$ distinct patches on proteins engaging in $`N/2`$ pairwise interactions of high specificity. These are not necessarily due to $`N`$ proteins. $`N`$ is an even number.

In practice, in a cell the pairwise binding interactions vary widely in strength, and some protein binding sites, which we model by patches, bind to more than one other protein. Our model protein can be used to deal with this variation but for simplicity we take all $`N/2`$ interactions to be equivalent. Also, as the number of proteins increases individual proteins become more dilute, this will tend to require larger equilibrium constants for binding. We neglect this effect here.

So, in a mixture with $`N`$ binding patches, each odd numbered patch $`i`$ is required to bind with the $`(i+1)`$th patch, with a large equilibrium constant $`K_b`$. Now, the $`i`$th ($`i`$ odd) and $`(i+1)`$th proteins must bind only to each other, it must be specific, and so in addition to binding to each other we require that the equilibrium constant for binding of patch $`i`$ to any patch other than patch $`i+1`$ be less than $`K_sK_b`$. Similarly, the equilibrium constant between patch $`i+1`$ and any protein other than protein $`i`$ must also be less than $`K_s`$. Thus if we denote the binding constant between proteins $`i`$ and $`j`$ by $`K_{ij}`$ we require

$`K_{ij}`$ $`>`$ $`K_b|ij|=1\mathrm{min}\{i,j\}\text{odd}`$
$`K_{ij}`$ $`<`$ $`K_s\text{otherwise}`$ (1)
## 3 Model protein

Our model protein is as first defined in , although there only solubility, not binding was considered, and because of this only patches with rather short sequences were considered. The model protein is chosen to be as simple and as generic as possible, while having interactions which are mediated by surface patches whose interactions are a function of the values of a sequence or string of elements. Figure 1 is a schematic of the model. Note that although the elements determine the interactions in our model and the residues determine the interactions in a real protein, there is no one-to-one relationship between one of our elements and an amino acid.

The model protein is a cube, with each of its 6 faces having a single patch, for details see . The lattice is cubic and each protein occupies 8 lattice sites arranged 2 by 2 by 2. We make the model 2 sites across to reduce the range of the attraction, which is 1 site, to half the diameter of the hard core: see  where an earlier version of this model was defined and its behaviour compared to that of real proteins. The model proteins can rotate, and so have 24 distinct orientations. Each of the 6 faces of the cube has a patch. If the faces of 2 proteins are in contact there is an energy of interaction between the 2 touching patches of the 2 proteins. By in contact we mean that the faces must overlap completely otherwise the energy of interaction is taken to be zero. Also, the model is such that the energy of interaction between two touching patches is a constant which does not change when the two proteins are rotated about the axis joining their centres.

How a patch interacts is specified by a sequence or string of $`n_E`$ elements. If a element has a value of 1 then the element is said to be hydrophilic or polar, whereas if it has a value of 0 then it is hydrophobic. A schematic of our model protein is shown in figure 1, where the sequence is represented by a ‘barcode’ with light stripes indicating 1s and dark stripes indicating 0s. The interaction energy, $`u_{ij}`$, of a pair of touching patches, $`i`$ and $`j`$, is then given by

$$u_{ij}=ϵ\underset{\alpha =1}{\overset{n_E}{}}\left(b_\alpha ^{(i)}1\right)\left(b_{1+n_E\alpha }^{(j)}1\right),$$
(2)
where $`b_\alpha ^{(i)}`$ is element number $`\alpha `$ of patch $`i`$. $`ϵ`$ is the interaction energy of 2 elements that are hydrophobic. We use energy units such that the thermal energy $`k_BT=1`$.

To calculate the interaction the string of elements of one of the patches is reversed and then the energy is just $`ϵ`$ times the sum of the number of pairs of corresponding elements where both elements are 0, are hydrophobic. The only interaction is between 2 hydrophobic elements; there is no hydrophobic-hydrophilic or hydrophilic-hydrophilic interaction. The reason one of the strings is reversed is that if this is not done then the interaction between like patches, $`j=i`$, is just $`ϵ`$ times the number of 0s in the string. Then all patches would stick to themselves. Reversing the strings removes this problem in a simple way. Of course, the interactions form a symmetric square matrix, $`u_{ij}=u_{ji}`$. Thus, the binding site of a protein is specified by giving values to the string of $`n_E`$ elements, and so there are $`2^{n_E}`$ possible different patches.

The binding constant $`K_{ij}`$ of a patch is given by

$$K_{ij}=\left(1/2\right)^{\delta _{ij}}\mathrm{exp}\left(u_{ij}\right)/6,$$
(3)
where the first factor is a symmetry factor which halves the equilibrium constant if $`i=j`$. The factor of $`1/6`$ is a normalisation factor of $`1/24`$, from the 24 distinct orientations, times the 4 possible rotations about the axis joining the centres of the proteins, all of which allow binding.

So, our interaction is mediated by patches, of which we have $`N`$ patches interacting in pairs. We neglect the other patches on the surface of the proteins, simply assuming that they are highly hydrophilic, having few hydrophobic elements, and so can be neglected. We do not need to specify the number of proteins $`N_{PR}`$, although of course as each protein has 6 patches, we must have that $`N_{PR}N/6`$.

Although we are focusing on proteins, a mixture of particles that satisfies equation (1) will self-assemble in a controlled way, the particles which have patch 1 will only stick to those that have patch 2, those with patch 3 stick to those with patch 4 and so on. We can even consider macroscopic objects with interaction energies determined by sequences made of large hydrophobic and hydrophilic stripes, e.g., a cube just as in figure 1 and perhaps a few mm or cm across. If these objects in water are agitated sufficiently strongly to break the weak bonds, those with $`K<K_s`$, but not the stronger bonds, those with $`K>K_b`$, then the cubes should self assemble. For systems roughly equivalent to our model with small $`n_E`$, this has been done in a series of beautiful experiments by Whitesides and coworkers .

## 4 Calculation of the capability

As each patch must recognise a unique other patch, each patch must have a unique sequence of hydrophobic and hydrophilic elements. As with $`n_E`$ elements there are $`2^{n_E}`$ different sequences then immediately we know that $`2^{n_E}N`$, or $`n_E\mathrm{log}_2(N)`$. This is an obvious lower bound on the number of elements required for proteins to bind specifically to 1 other protein from $`N`$. It does not take into account our requirement, equation (1). If we use all $`2^{n_E}`$ sequences then many sequences in the mixture will differ from other sequences present by the value of only 1 element. Yet we require that a patch on a protein bind very strongly to one other patch and very weakly to all the others. For example, we require that patch 1 bind to patch 2 strongly but patch 3 very weakly or not at all. Thus, we require that $`K_{12}>K_b`$ and $`K_{13}<K_s`$. But if say proteins 2 and 3 differ in only one element then $`K_{12}`$ and $`K_{13}`$ can differ by a factor of at most $`\mathrm{exp}(ϵ)`$. So, then $`K_{12}/K_{13}\mathrm{exp}(ϵ)`$, which is inconsistent with our requirements for $`K_{12}`$ and $`K_{13}`$ unless $`ϵ`$ is large. In other words as the ratio $`K_b/K_s`$ is large, in order satisfy our criterion for functionality, equation (1), we cannot have 2 sequences present in which only 1 element is different. This requirement to ‘space out’ the patches dramatically reduces the number of proteins we can simultaneously use in a mixture. We will need considerably more than $`\mathrm{log}_2(N)`$ elements per protein patch to achieve a functional mixture.

In fact this problem of a patch picking out one other patch from a total of $`N`$ patches is analogous to problems that arise in the study of communication through noisy channels. There the receiver of a message must pick out which of the possible messages they have received. The study of communication channels gave birth to what is now often called information theory . In the 1940s Shannon wanted to know how many messages could be sent through a given noisy channel during some time interval, with each message being received correctly . Reference  is reprinted in . The theory is statistical in nature, in the sense that it does not consider a specific message but considers typical messages. Our theory too is statistical, we will not consider a specific mixture.

Now, we require that a protein bind to another with a large equilibrium constant $`K_b`$. This will require that the protein have a minimum of $`B`$ hydrophobic elements, where

$$B=\text{int}\left(\mathrm{ln}(6K_b)/ϵ\right)+1,$$
(4)
from equation (3). The function $`\text{int}(x)`$ yields the largest integer less than $`x`$. It is required as the number of elements $`B`$ must of course be an integer. As we will see, the more hydrophobic elements a patch has the more other patches it sticks to, i.e., the more other possible patches have $`K_{ij}>K_s`$, and so the number of hydrophobic elements should be kept to a minimum. We demonstrate this at the end of this section. As it cannot be less than $`B`$ while still binding strongly to one other patch, the optimum number of hydrophobic elements is $`B`$. From now on, we only consider patches with $`B`$ hydrophobic elements. With this restriction the number of possible patches, $`N_{POS}`$, is

$$N_{POS}=\begin{array}{c}\text{Number of}\\ \text{proteins with }B\\ \text{hydrophobic elements}\end{array}=\frac{n_E!}{(n_EB)!B!}.$$
(5)
Actually, this neglects the fact that some of these sequences are not actually possible as they bind to themselves with $`K_{ii}>K_s`$, and this makes equation (5) a slight overestimate. The mean number of interactions of a patch with itself is close to $`B^2/(2n_E)`$. In the next section we will be considering patches where $`B^2/(2n_E)`$ is fewer than the number of hydrophobic interactions that make an interaction too sticky. So we will use equation (5) and neglect the fact that a small fraction of patches stick to themselves.

But any one of the possible proteins will stick to many other proteins. We use the term ‘stick to’ to mean bind with an equilibrium constant $`K_{ij}>K_s`$. We need to calculate the number of proteins $`j`$ for which $`K_{ij}>K_s`$. We denote this number by $`N_S`$. If $`S`$ is the minimum number of elements for which $`K_{ij}>K_s`$, then $`S`$ is given by

$$S=\text{int}\left(\mathrm{ln}(6K_s)/ϵ\right)+1.$$
(6)
Now, if $`N_s`$ is the number of patches with $`B`$ hydrophobic elements, $`s`$ of which interact with a fixed set of $`B`$ elements on another patch, then the number of patches that stick to a given patch, $`N_S`$, is

$$N_S=\underset{s=S}{\overset{B}{}}N_s,$$
(7)
because any patch with $`sS`$ hydrophobic interactions contributes to the number of patches that are too sticky. Calculating $`N_s`$ is a simple exercise in combinatorics. Consider a patch with $`B`$ elements. The number of ways, $`N_s`$, that $`s`$ elements can be situated such that they interact with a fixed pattern of $`B`$ elements out of $`n_E`$ is

$$N_s=\frac{B!}{(Bs)!s!}\frac{(n_EB)!}{(n_EB(Bs))!(Bs)!}.$$
(8)
Having found expressions for the total number of patches, $`N_{POS}`$, and the number each one sticks to, $`N_S`$, we proceed to see how many specific interactions are possible. The first thing to note is that once a patch is specified, so is its partner. All $`B`$ of the partner’s hydrophobic elements must be at the $`B`$ positions that interact with the hydrophobic elements of the first protein and so there is no freedom in choosing its sequence of elements. Thus, as there is a unique partner for each patch, we have $`N_{POS}/2`$ binding pairs of patches. Actually a very small fraction of patches bind to themselves, but we neglect this.

So, we have $`N_{POS}/2`$ pairs. Every time we add one of these pairs to the mixture, we have to discard approximately $`2N_S`$ of these pairs. This figure of $`2N_S`$ is obtained by making two assumptions. The first is that each of the 2 patches of the pair eliminates $`N_S`$ pairs by sticking to one of the partners — we neglect sticking to more than one patch of a binding pair. The second is that there is no overlap between the sets of approximately $`N_S`$ pairs that each of the 2 patches of the pair added sticks to. If every time we add a pair we have $`2N_S`$ fewer pairs then clearly we run out of pairs after we have added $`(N_{POS}/2)/(2N_S)`$ pairs. We name the maximum number of patches that can bind together in pairs in a mixture, the capability and denote it by $`C`$. For our model it is

$$C(n_E,ϵ)=\frac{N_{POS}}{2N_S}.$$
(9)
The maximum number of high specificity bonds is of course half the capability $`C`$. In information theory, the fundamental quantity is the capacity $`C_{Sh}`$. In units of bits per bit sent, it is related to our capability by $`C_{Sh}=\mathrm{log}_2(C)/n_E`$ . In communication the analogue of our capability $`C`$ increases exponentially so it is convenient to use not the capability itself but its logarithm the capacity. As we will see our $`C`$ does not increase exponentially with $`n_E`$ and so it is more convenient to work with it directly. In an appendix we discuss the similarities and differences to communication.

To derive equation (9) we assumed that $`n`$ patches stick to $`nN_S`$ others, neglecting any overlap between the $`N_S`$ patches that stick to one patch and the $`N_S`$ patches that stick to another. A simplifying assumption that means that equation (9) will tend to underestimate the true capability.

For their pioneering study of the immune system, Perelson and Oster  introduced the idea of ‘shape space’ of proteins where if a protein from the immune system is near in this space to a foreign protein, it will bind to, recognise, the foreign protein. Our model proteins too have a corresponding space and the capability we have calculated essentially corresponds to the number of proteins required to fill this space. By fill we mean have so many proteins that any extra one of the $`N_{POS}`$ patches with $`B`$ hydrophobic elements that is added is likely to stick to one of the patches in the mixture. For an immune system to function the space should be overfilled to avoid any foreign protein escaping detection , it is an advantage for more than one immune-system protein to bind to a foreign protein, whereas for specific binding the space must be underfilled.

When the first term in equation (7) for $`N_S`$ is the dominant one, we can approximate $`N_S`$ as

$`N_S`$ $``$ $`{\displaystyle \frac{B!}{(BS)!S!}}{\displaystyle \frac{(n_EB)!}{(n_EB(BS))!(BS)!}}`$ (10)
$``$ $`{\displaystyle \frac{B!n_E^{BS}}{S!(BS)!^2}},`$
where in the second line we simplified the expression assuming $`n_EB`$ is not too small. Making similar approximations for $`N_{POS}`$ we get $`N_{POS}n_E^B/B!`$ and hence for the capability $`C`$,

$`C`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{n_E^B}{B!}}{\displaystyle \frac{S!(BS)!^2}{B!n_E^{BS}}}={\displaystyle \frac{(BS)!^2S!}{2B!^2}}n_E^S.`$ (11)
Also, note that from this approximate expression it is clear that to maximise the capability the number of hydrophobic elements, $`B`$ in equation (11) should be as small as possible. In equation (11), the $`B`$ dependence is $`C1/[B(B1)\mathrm{}(BS+1)]`$ and this is a monotonically decreasing function of increasing $`B`$. Finally, as promised $`C`$ does not increase exponentially with the number of elements $`n_E`$. We should note that Shannon’s channel-capacity theorem  implies that there potentially exists a model in which $`C`$ does increase exponentially. So real proteins may have a $`C`$ that increases exponentially with the number of amino acids in a binding patch, however in our model the increase is sub-exponential.

## 5 Results for the capability

The binding strength of protein interactions are often given in terms of the dissociation constant, typically with Molar (M) units . 1 Molar is $`6\times 10^{23}`$ molecules per litre, which is $`0.6`$ molecules/nm<sup>3</sup>. The dissociation constant, $`K_d`$, for $`A`$ and $`B`$ binding together is, by definition, $`K_d=c_Ac_B/c_{AB}`$, where $`c_A`$, $`c_B`$ and $`c_{AB}`$ are the Molar concentrations of free (unbound) $`A`$, free $`B`$ and bound $`AB`$ dimers, respectively . It is therefore just one over the equilibrium constant for binding. Dissociation constants of protein interactions vary over orders of magnitude. They can easily be of order $`\mu `$M and in some cases values as low as $`10^{16}`$M have been measured. For example, some E. coli can produce a protein called E9 (it is one of a family of related proteins) which is highly toxic, it is an antibiotic. E9 is neutralised by another protein, Im9, which binds to it to form a heterodimer . At low salt in vitro the dissociation constant for the binding of Im9 to E9 is close to $`10^{16}`$M. In $`250`$mM sodium chloride solution the dissociation constant is a mere $`10^{14}`$M: the binding has a significant electrostatic component. This is a remarkably strong and specific binding, it must be specific as it must avoid binding in vivo to proteins other than E9.

E9 and Im9 are a somewhat extreme example, to get a feel for more typical values, consider a bacterium such as E. coli. It has a cytoplasm with a volume around 1$`\mu `$m<sup>3</sup>. So, one protein in this volume has a concentration of $`10^9`$nm<sup>-3</sup> or about 1nM. If there are only say 10 copies of $`A`$ and $`B`$ then for a significant fraction of these copies to bind to each other, then a dissociation constant around a $`10`$nM is required. If there are around 1000 copies of the proteins $`A`$ and $`B`$ in each E. coli cell then a dissociation constant of around 1$`\mu `$M is enough to ensure roughly equal numbers of the free $`A`$ and $`B`$ and of $`A`$ and $`B`$ bound into dimers.

Our proteins are 2 lattice sites across and our equilibrium constants $`K_{ij}`$ are in units of the volume of a lattice site. Proteins are a few nms across and so each lattice site corresponds to a volume about 1 to 2nms across. Taking a lattice site to be 1nm across, and noting that 1Molar is close to 1 per nm<sup>3</sup> we have that $`K_{ij}=1`$ is approximately equivalent to an equilibrium constant of 1M<sup>-1</sup>. Thus to obtain a dissociation constant of 1$`\mu `$M we need a $`K_{ij}10^6`$, and so need to set $`K_b=10^6`$.

In figure 2 we have plotted the capability $`C`$ for $`K_b=10^6`$ (solid curve) and $`10^8`$ (dashed curve). In both cases we insisted that the equilibrium constant between patches not binding to each other be less than $`K_s=100`$. A $`K_b=10^8`$ is approximately equivalent to a dissociation constant of 10nM. We took the energy of interaction between hydrophobic elements $`ϵ=1`$. Then for $`K_b=10^6`$, $`B=16`$ and, for $`K_b=10^8`$, $`B=21`$, and with $`K_s=100`$, $`S=7`$ elements is too sticky, $`K_{ij}>100`$. Also, the dotted curves are the result of equation (11).

Of course the capability $`C`$ increases as $`n_E`$ increases but for the stronger binding, even for $`n_E=150`$, $`C=40`$. This is highly inefficient, the restriction to be functional means that we need $`150`$ elements for a mixture with 40 specific binding interactions, whereas if we just needed 40 different patches then $`\mathrm{log}_240=5.3`$ elements would have been enough. However, each patch of a protein is required to bind to 1 other patch with a very large equilibrium constant, and with equilibrium constants at least 6 orders of magnitude smaller to any of the other 39 patches. This requirement is a very stringent one.

The inefficiency can be clearly seen in figure 3 where we have plotted the number of possible patches to which a patch sticks. We see that if we insist on specific binding of strength $`K_b=10^8`$ (the dashed curve) then we require $`B=21`$ hydrophobic elements for bonding of this strength but then when $`n_E=150`$, a patch with 21 hydrophobic elements sticks to approximately $`2.8\times 10^{23}`$ other possible proteins. The equilibrium constant between a single patch and of order $`10^{23}`$ others is greater than $`K_s=100`$. Thus to make a mixture that is functional according to our criteria, for every specific-binding patch we add to the mixture the number of possible patches that we can no longer add to the mixture increases by $`10^{23}`$.

The same thing must apply to real proteins in vivo, although proteins are much more complex than our, deliberately very simple, model and so it is not possible to quantify $`N_S`$ for real proteins. However, we can return to our example of E9 and Im9 . A number of single-residue mutants have been made for approximately 30 of the residues of Im9 suspected of contributing to the E9-Im9 binding. The changes in the dissociation constant in each case were measured, and they cover a large range  but many of them resulted in only modest increases in the dissociation constant. Thus, many sequences with these mutations are not viable as part of the amino-acid sequences of proteins. We cannot have, for example, an enzyme that has a part of its amino-acid sequence close to that of Im9 as the enzyme will bind to E9. Many more mutants, mutants with 2, 3, etc. mutations, will also bind strongly to E9 and in general this will be inappropriate. So functional proteins cannot include these amino-acid sequences, thus reducing the number of amino-acid sequences which are functional in the mixture of proteins in the cytoplasm of E. coli when E9 is present.

## 6 Conclusion

The mixture of proteins inside cells is not simply an inert mixture like, for example, the complex mixture of hydrocarbons in crude oil. It has evolved to be functional: the proteins do things. In particular, many proteins must bind strongly to a specific patch or patches on one or more other proteins. Here, we studied a very simple model protein, introduced in . We borrowed ideas from the theory of communication, and calculated the maximum number of specific binding interactions that can function together in a mixture. We called it the capability $`C`$. It is roughly the number of possible different patches, $`N_{POS}`$, over the number of patches that stick to a given patch and so cannot be in the same mixture, $`N_S`$, $`C=N_{POS}/2N_S`$, The capability $`C`$ of the constituent proteins limits the size of the interactome. The interactome is a network or web of many specific interactions, and it can have at most $`C/2`$ links. A network with many links requires that the proteins of which it is composed must have large patches, patches with a large number of elements.

Schneider has applied information theory to analysing the DNA sequences of sequences of base pairs to which proteins bind with high specificity . We briefly apply an analogous analysis to our model in the appendix. Work on the size of the regions on foreign proteins that are recognised by the immune-system proteins is very similar in approach to that used here . For example, Percus et al.’s  study of how demanding is the requirement that our immune-system proteins differentiate between foreign and human proteins, is analogous to our study of how demanding is the requirement that a protein bind to only one other patch in a mixture. However, we do use a concrete although highly simplified, model of a protein unlike work on the immune system which has used rather abstract models of recognition. Indeed, it might be useful to use our concrete model of a protein in future studies of the immune system.

The more specific the binding, the larger $`K_b/K_s`$ at fixed $`K_s`$, the more elements $`n_E`$ we need to achieve a given capability $`C`$, see figure 2. Now, although there is no one-to-one correspondence between our hydrophobic/hydrophilic elements and the residues of a real protein, we expect that increasing $`n_E`$ corresponds in real proteins to increasing the size of the protein. But ideally, proteins should be as small as possible, they will then consume less material (amino acids), diffuse faster, and take up less volume. For these reasons we might want $`n_E`$ to be as small as possible and so there is a trade off between the size and complexity needed for specific binding and the problems of being large. This trade off involves not only the properties of an individual protein but also the global properties of the mixture, e.g., the number of components. Eukaryote cells are compartmentalised  and each compartment will have only a subset of the total number of proteins. A smaller number of proteins with higher concentrations in a compartment allows strong specific binding at smaller values of $`n_E`$.

Very recent work on the specificity of protein-protein binding has been done by Havranek and Harbury , and by Shifman and Mayo . They looked at specificity but in both cases they used a combination of a much more detailed model and experiments. In particular, Havranek and Harbury considered 2 proteins and constrained each protein to: i) to dimerise with itself, ii) not to bind to the other protein, iii) not to aggregate. They found that proteins have to be specifically engineered not to bind to other proteins and not to aggregate, as well as to bind to their binding partner (themself in this case). The detailed nature of their model restricted them to a mixture of 2 proteins. It also precluded the exhaustive enumeration of all possible proteins and elimination of all of the proteins that bind to proteins other than their partners, that we performed here. Following the approach taken by Shannon  we also did not consider a specific mixture, as Havranek and Harbury did, but calculated a limit on what the mixture could do, here the limit on the number of interaction of a given specificity. However, theirs and our approaches are related. We both consider how specificity is achieved and in doing so explicitly consider the binding of a protein to the wrong other protein or site, that must be prevented. This avoidance of an undesired state is referred to as negative design .

Future work could also go in a number of other directions. The concepts introduced here of a capability $`C`$ for specific binding of a protein could be applied to other models of proteins. Also, the proteins in a cell are the product of evolution of course, and it would be interesting to study the constraints on evolution of proteins due to the requirement of specific binding. Very recent work  on large data sets of interactions, partial interactomes, has found a negative correlation between the rate of evolution of the amino-acid sequence of a protein and the number of specific interactions it participates in. Although there is some debate over this finding . This work relies on recent work characterising very large numbers of interactions, particularly in yeast . Also see the references of  for earlier work on individual proteins. Within our model, if a protein patch is part of a mixture near the capability of its proteins, then not only would many mutations be deleterious because they weakened its binding to its partner, many mutations would be deleterious because they caused the patch to stick to other patches in the mixture. The mutation rate of an individual protein will in general depend on the sequences of elements of all the other patches, not just that of its partner.

It is a pleasure to acknowledge that this work started with inspiring discussions with D. Frenkel. I would also like to acknowledge discussions with J. Cuesta, and J. Doye and A. Louis for introducing me to negative design. Finally, I would like to thank an anonymous reviewer for drawing my attention to the connection to the work of Whitesides. This work was supported by the Wellcome Trust (069242).

## Appendix

In our notation, the channel-capacity theorem of Shannon states that in the $`n_E\mathrm{}`$ limit, we can achieve a capability $`C`$ such that $`\mathrm{log}_2(C)/n_EO(1)`$, with fixed $`K_s`$ and the specificity ratio $`K_b/K_s`$ arbitrarily large. For our model protein $`C`$ does not increase exponentially with $`n_E`$ and so in the $`n_E\mathrm{}`$ limit, $`\mathrm{log}_2(C)/n_E0`$. Thus the channel-capacity theorem does not apply to our model of a protein. Note that this does not rule out a different type of model protein for which the channel capacity theorem holds, or that real proteins might obey it.

Shannon’s derivation of the channel-capacity theorem (theorem 11 of ) breaks down when applied to our proteins because he assumes, correctly for communication, that the entropy per bit sent of the strings is some constant which he calls $`H`$. Here the entropy per element $`s`$, is

$$s=\left(\frac{B}{n_E}\right)\mathrm{ln}\left(\frac{B}{n_E}\right)\left(\frac{n_EB}{n_E}\right)\mathrm{ln}\left(\frac{n_EB}{n_E}\right)$$

which is the logarithm of equation (5) over $`n_E`$. At fixed $`ϵ`$, $`B`$ is fixed and so as $`n_E`$ is increased, $`s`$ decreases. In order to avoid the patches becoming very sticky the number of hydrophobic elements must be minimised and this makes $`s`$ decrease with increasing $`n_E`$ which in turn dramatically reduces the capability. This feature is present in our model but not in communication systems. It is also not present in the more abstract approach to applying information theory to high specificity binding of Schneider . It is possible that it is simply an artifact of our simple model, but other related models do show a similar behaviour . Another perspective on this is to note that to make $`s`$ not a function of $`n_E`$, $`B`$ should not be a constant, it should be $`fn_E`$, with $`f<1`$ a constant. But then the average interaction energy between 2 patches is $`ϵf^2n_E/2`$: which increases with $`n_E`$. In essence the problem is that the larger the number of elements of the patch the larger the interaction energy is. It is possible that a similar problem also occurs in real proteins, where if a large patch on the surface of a protein is required to achieve binding of high affinity it is difficult to avoid this large patch interacting strongly with many other proteins. The analogue of the $`\mathrm{log}(N_S)/n_E`$ in communication, called the equivocation, is also a constant, unlike $`\mathrm{log}(N_S)/n_E`$.

Finally, it should be noted that $`sn_E`$ is simply the entropy associated with the space of proteins that have the right number of hydrophobic elements. This is the space consisting of $`N_{POS}`$ proteins. It is not an information in the sense that the information of binding sites on DNA evaluated by Schneider  is an information. He evaluates the information about a DNA sequence that you obtain if you are told that the sequence is the binding site for say the protein LacI . This information, call it $`I_s`$, is, by definition, the difference between the entropy of the DNA sequence ($`2\times \text{number of base pairs}`$ if the 4 base pairs all occur with probability $`0.25`$) if you know nothing about it and the entropy of sequences that LacI binds too. We can evaluate a similar information, again call it $`I_s`$, for our model proteins. Let us consider protein $`i`$, with $`i`$ odd so that it binds to protein $`i+1`$. $`I_s`$ is the information about the sequence of protein $`i`$ that you obtain if you are told what the sequence is of protein $`i+1`$. In section 4 when we calculated the capability $`C`$ we insisted on proteins having $`B`$ hydrophobic elements and then once the sequence of protein $`i+1`$ is fixed, so is that of protein $`i`$ and so the entropy of sequence $`i`$ is then 0 once the sequence of $`i+1`$ is known. Then $`I_s=sn_E`$. But $`I_s`$ need not always equal $`sn_E`$. In a mixture of proteins with interactions coded by $`n_E`$ elements, with $`N`$ below $`C(n_E)`$ more than one sequence $`i`$ will bind to $`i+1`$ and also satisfy the requirement not to stick to the other $`N1`$ proteins. The number of hydrophobic elements will not always be $`B`$. Then the entropy of protein $`i`$ once the sequence of protein $`i+1`$ has been specified will be greater than 0 and $`I_s`$ will not equal the entropy of the space of possible proteins.

## References


text: # Withdrawn


text:
# Single Crystal Growth of Skutterudite CoP3 under High Pressure

## 1 Introduction

Filled skutterudite compounds, RM4X12 (R = rare earth, M = Fe, Ru or Os; X = P, As or Sb), are of interest to the scientific community because of their potential as thermoelectric materials and the instability of their $`f`$ electrons. It is believed that the $`f`$ electron instability contributes to a wide variety of phenomena in the skutterudites, such as superconductivity, a metal-insulator transition  and magnetic ordering.

To clarify the nature of the $`f`$ electron instability, as well as to develop the thermoelectric performance of skutterudites, studies based on single crystals are essential. In particular, to improve the thermoelectric properties, it is important to clarify the origin of their superior thermoelectric performance. Filled skutterudites exhibit a remarkably low thermal conductivity, which enhances their performance. Great efforts, thus, have been made to clarify the origin of this low thermal conductivity. One proposed model that could explain the phenomenon is the so-called rattling effect. It is conjectured that phonons are scattered effectively by the free vibration of rare earth atoms in large lattice cages . To examine the rattling effect, phonon behavior should be studied by neutron scattering, which requires large single crystals.

To date, a number of groups have grown single crystals of various skutterudites. According to their reports, skutterudites are incongruent melting compounds. In the growth process, Sb or Sn were typically used as flux, and samples were sealed in evacuated silica tubes. The volume of grown binary skutterudites, such as CoSb3 was typically 0.6 cm3 . For ternary filled skutterudites, on the other hand, a typical sample size was about $`1\times 1\times 1`$ mm3 , considerably smaller than that of binary skutterudites.

In the growth of skutterudite single crystals by the flux method, however, several limitations have been identified. First, only a few skutterudite compounds have been able to be crystallized. This is a disadvantage in a comprehensive study of physical properties, where single crystals of a number of different skutterudites are necessary. Second, the size of the grown crystals of ternary filled skutterudites is limited. Larger single crystals are necessary in order to study fundamental properties by some methods, for example, by neutron scattering measurements. To solve these problems, the development of an original approach to the growth of skutterudite crystals may be beneficial.

In this study, we analyze the growth of single crystals of skutterudites performed under high pressure using a wedge-type, cubic-anvil, high-pressure apparatus. This technique was previously applied to enhance the growth of black phosphorus single crystals . Using this method, single crystals of black phosphorus of dimension $`4\times 2\times 0.2`$ mm3 were successfully grown under a pressure of 2.3 GPa. In this paper, we report that this technique is also useful for growing single crystals of skutterudites.

## 2 Experimental Details

CoP3 single crystals were grown using a wedge-type, cubic-anvil, high-pressure apparatus (CAP-07, RIKEN). Fig. 1 shows a schematic representation of the sample cell assembly used in the experiment. The sample container is a cube composed of pyrophyllite with sides of length 15.3 mm. A carbon heater tube and a BN (boron nitride) crucible were inserted into the container. The size of the BN crucible was 6 mm in length and 4.5 mm in diameter. To estimate the temperature in the crucible during crystal growth, the relationship between the applied electrical power and the temperature was predetermined without encasing the samples. Temperatures were measured up to $`1100^{}`$C at the center of the carbon tube using Chromel-Alumel thermocouples at ambient pressure without employing any emf correction for the pressure effect. Temperatures above $`1100^{}`$C were estimated by extrapolating the obtained temperature calibration curve.

The raw materials of Co and P of 99.9crucibles after mixing them in a stoichiometric ratio. The filled crucibles were then immediately placed in the container and compressed to 3.5 GPa at room temperature. The samples were then heated to $`1100^{}`$C and sintered for 1 hour to pre-synthesize CoP3 polycrystals. Subsequently, the samples were heated to $`1500^{}`$C and maintained at that temperature for 1 hour. Finally, the samples were cooled to $`1000^{}`$C over a period of 5 hours, followed by quenching to room temperature before pressure release.

The grown crystals were characterized by the following methods. Phase identification of the grown crystals was performed by powder x-ray diffraction (RINT-1000, RIGAKU) at room temperature using Cu Ka radiation. The composition of the grown crystals, as well as the distribution of the Co and P atoms, was examined by a JEOL JSM-6301F scanning electron microscope x-ray energy dispersion spectrometer (SEM-EDS) with an accelerating voltage of 20 kV. A Li doped Si crystal was used as the SEM-EDS detector with an ultra-thin polymer window (CDU-Super-UTW, EDAX), that could detect elements ranging from .Be to U. To estimate the concentration of Co and P atoms, characteristic x-rays of Co Ka and P Ka were used. A mapping analysis of the Co and P atoms was conducted with a scanning time of 10 ms/point and a measurement point density of 2500 points/mm2 over the entire boule cross-section. For the SEM-EDS measurements, the as-grown crystals were mechanically polished to mirror-like-surfaces. Back-reflection Laue x-ray measurements were performed with a collimator 1 mm in diameter to confirm that the grown crystals were single crystals.

## 3 Results

Fig. 2(a) shows a cross section of the as-grown boule of CoP3. The shape of the boule is cylindrical, about 3.4 mm in diameter, with a number of cracks. The as-grown boule split quite easily along the radial direction to generate specimens $`1`$ mm in diameter and $`0.3`$ mm in length (Fig. 2(b)). Fig. 3 shows a back-reflection Laue x-ray photograph of the specimen displayed in Fig. 2(b). In the photograph, 2mm symmetry is observed, indicating a (100) crystallographic plane. This result suggests that the specimen is a single crystal.

Powder x-ray diffraction measurements were carried out using powder ground from the as-grown boule (Fig. 4). The observed diffraction lines are all

indexable using the skutterudite structure (space group Im$`\overline{3}`$). No impurity phases are observed, suggesting that the entire as-grown boule is a single phase of CoP3. The lattice constant determined by a least squares fit to the data was a = 7.7067 , which is close to the literature value . This agreement suggests that the grown crystals are of high purity.

Fig. 5 shows the distribution of Co and P atoms within the as-grown boule measured by SEM-EDS. The presence of the darker regions in the map is due to cracks, presumably formed during the cooling process. In the other regions, a homogeneous distribution of Co and P atoms is observed. We checked carefully that there were no additional elements present, especially around the cracks. The ratio of the Co to the P atoms was determined by averaging several measurements and was found to be Co:P = 0.97:3.00, which is consistent with a stoichiometric composition within instrumental error. These results suggest that the crystals are a high purity, single phase of CoP3.

## 4 Discussion

In this study, single crystals of CoP3 were grown from a stoichiometric melt under high pressure using a wedge-type, cubic-anvil, high-pressure apparatus. The entire as-grown boule was characterized as a single phase of CoP3 by powder x-ray diffraction and SEM-EDS measurements. These results suggest that CoP3 is a congruent melting compound under high pressure. Around ambient pressures, however, it is still unclear whether CoP3 is a congruent melting compound or not, since the phase diagram of Co-P is incomplete . Note that other skutterudites whose phase diagrams are well known, for example CoAs3 and CoSb3, are all incongruent melting compounds around ambient pressures . It is, thus, natural to assume that CoP3 is also an incongruent melting compound around ambient pressures in contrast to the results under high pressure. Most likely, CoP3 transforms from an incongruent to a congruent melting compound by the application of high pressures. If this presumption is true, other skutterudites may also melt congruently under high pressure, which would be a great advantage for growing single crystals.

The size of the CoP3 single crystals produced in this study is typically $``$ 1 mm in diameter and $`0.3`$ mm in length. To enlarge the size of single crystals, the cracks in the as-grown boules need to be eliminated. One possible solution to prevent cracking is to use a slower cooling rate. Another solution could be to apply a more homogeneous pressure by using a smaller crucible. Crystal growth from a single seed is also important for obtaining large single crystals. For this to take place, the temperature distribution in the furnace must be further optimized.

## 5 Conclusion

Single crystals of CoP3 were grown from a stoichiometric melt under a pressure of 3.5 GPa using a wedge-type, cubic-anvil, high-pressure apparatus. The size of the grown single crystals is typically $`1`$ mm in diameter and $`0.3`$ mm in length. The as-grown boules were characterized by powder x-ray diffraction and SEM-EDS measurements and were determined to be a single phase of CoP3. These results suggest that CoP3 melts congruently under high pressure.

## 6 Acknowledgements

The authors thank A. Iyo for valuable technical suggestions. We also thank A. Negishi and A. Yamamoto for help in performing SEM-EDS measurements. This work was supported by a Grant from the Ministry of Economy, Trade and Industry of Japan.