text:
# Mapping class groups are linear

## 1 Introduction

A. The Harvey conjecture. The mapping class group has been introduced in the 1920-ies by M. Dehn . Such a group, $`Mod(X)`$, is defined as the group of isotopy classes of the orientation-preserving diffeomorphisms of a two-sided closed surface $`X`$ of genus $`g1`$. The group is known to be prominent in algebraic geometry , topology  and dynamics . When $`X`$ is a torus, the $`Mod(X)`$ is isomorphic to the group $`SL_2(Z)`$. (The $`SL_2(Z)`$ is called a modular group, hence our notation for the mapping class group.) A little is known about the representations of $`Mod(X)`$ beyond the case $`g=1`$. Recall, that the group is called linear, if there exists a faithful representation into the matrix group $`GL_m(R)`$, where $`R`$ is a commutative ring. The braid groups are known to be linear . Using a modification of the argument for the braid groups, it is possible to prove, that $`Mod(X)`$ is linear in the case $`g=2`$ . Whether the mapping class group is linear for $`g3`$, is an open problem, known as the Harvey conjecture , p.267.

B. The Teichmueller functor. A covariant (non-injective) functor from a category of generic Riemann surfaces to a category of the so-called toric $`AF`$-algebras (see section 2 for a definition) was constructed in . The functor, a Teichmueller functor, maps any pair of isomorphic (i.e. conformal) Riemann surfaces to a pair of the stably isomorphic (Morita equivalent) toric $`AF`$-algebras. Since each isomorphism of Riemann surfaces is given by an element of $`Mod(X)`$ , it is natural to ask about a representation of $`Mod(X)`$ by the stable isomorphisms of toric $`AF`$-algebras; roughly, our objective can be stated as follows.

Main problem. To study the Harvey conjecture from the standpoint of toric $`AF`$-algebras.

The stable isomorphisms of toric $`AF`$-algebras are well understood and surprisingly simple; provided the automorphism group of the algebra is trivial (this is true for a generic algebra), its group of stable isomorphism admits a faithful representation into a matrix group over the commutative ring $`Z`$ . This fact, combined with the properties of the Teichmueller functor, implies an amazingly simple positive solution to the Harvey conjecture.

###### Theorem 1

For every surface $`X`$ of genus $`g2`$, there exists a faithful representation $`\rho :Mod(X)GL_{6g6}(Z)`$.

The structure of the note is as follows. In section 2, the preliminary facts, necessary to prove theorem 1, are brought together. Theorem 1 is proved in section 3.

## 2 Preliminaries

We review the toric $`AF`$-algebras and the Teichmueller functor on the space of generic Riemann surfaces; the reader is referred to , ,  and  for details.

### 2.1 $`AF`$-algebras

A. The $`C^{}`$-algebras. By the $`C^{}`$-algebra one understands a noncommutative Banach algebra with an involution. Namely, a $`C^{}`$-algebra $`A`$ is an algebra over the complex numbers $`C`$ with a norm $`aa`$ and an involution $`aa^{},aA`$, such that $`A`$ is complete with the respect to the norm, and such that $`abab`$ and $`a^{}a=a^2`$ for every $`a,bA`$. If $`A`$ is commutative, then the Gelfand theorem says that $`A`$ is isometrically $``$-isomorphic to the $`C^{}`$-algebra $`C_0(X)`$ of the continuous complex-valued functions on a locally compact Hausdorff space $`X`$. For otherwise, the algebra $`A`$ represents a noncommutative topological space.

B. The stable isomorphisms of $`C^{}`$-algebras. Let $`A`$ be a $`C^{}`$-algebra deemed as a noncommutative topological space. One can ask, when two such topological spaces $`A,A^{}`$ are homeomorphic? To answer the question, let us recall the topological $`K`$-theory. If $`X`$ is a (commutative) topological space, denote by $`V_C(X)`$ an abelian monoid consisting of the isomorphism classes of the complex vector bundles over $`X`$ endowed with the Whitney sum. The abelian monoid $`V_C(X)`$ can be made to an abelian group, $`K(X)`$, using the Grothendieck completion. The covariant functor $`F:XK(X)`$ is known to map the homeomorphic topological spaces $`X,X^{}`$ to the isomorphic abelian groups $`K(X),K(X^{})`$. Let $`A,A^{}`$ be the $`C^{}`$-algebras. If one wishes to define a homeomorphism between the noncommutative topological spaces $`A`$ and $`A^{}`$, it will suffice to define an isomorphism between the abelian monoids $`V_C(A)`$ and $`V_C(A^{})`$ as suggested by the topological $`K`$-theory. The role of the complex vector bundle of the degree $`n`$ over the $`C^{}`$-algebra $`A`$ is played by a $`C^{}`$-algebra $`M_n(A)=AM_n`$, i.e. the matrix algebra with the entries in $`A`$. The abelian monoid $`V_C(A)=_{n=1}^{\mathrm{}}M_n(A)`$ replaces the monoid $`V_C(X)`$ of the topological $`K`$-theory. Therefore, the noncommutative topological spaces $`A,A^{}`$ are homeomorphic, if the abelian monoids $`V_C(A)V_C(A^{})`$ are isomorphic. The latter equivalence is called a stable isomorphism of the $`C^{}`$-algebras $`A`$ and $`A^{}`$ and is formally written as $`AKA^{}K`$, where $`K=_{n=1}^{\mathrm{}}M_n`$ is the $`C^{}`$-algebra of compact operators. Roughly speaking, the stable isomorphism between the $`C^{}`$-algebras means that they are homeomorphic as the noncommutative topological spaces.

C. The $`AF`$-algebras. An $`AF`$-algebra (approximately finite $`C^{}`$-algebra) is defined to be the norm closure of an ascending sequence of the finite dimensional $`C^{}`$-algebras $`M_n`$’s, where $`M_n`$ is the $`C^{}`$-algebra of the $`n\times n`$ matrices with the entries in $`C`$. Here the index $`n=(n_1,\mathrm{},n_k)`$ represents a semi-simple matrix algebra $`M_n=M_{n_1}\mathrm{}M_{n_k}`$. The ascending sequence mentioned above can be written as $`M_1\stackrel{\phi _1}{}\mathrm{M}_2\stackrel{\phi _2}{}\mathrm{},`$ where $`M_i`$ are the finite dimensional $`C^{}`$-algebras and $`\phi _i`$ the homomorphisms between such algebras. The set-theoretic limit $`A=limM_n`$ has a natural algebraic structure given by the formula $`a_m+b_ka+b`$; here $`a_ma,b_kb`$ for the sequences $`a_mM_m,b_kM_k`$. The homomorphisms $`\phi _i`$ can be arranged into a graph as follows. Let $`M_i=M_{i_1}\mathrm{}M_{i_k}`$ and $`M_i^{}=M_{i_1^{}}\mathrm{}M_{i_k^{}}`$ be the semi-simple $`C^{}`$-algebras and $`\phi _i:M_iM_i^{}`$ the homomorphism. One has the two sets of vertices $`V_{i_1},\mathrm{},V_{i_k}`$ and $`V_{i_1^{}},\mathrm{},V_{i_k^{}}`$ joined by the $`a_{rs}`$ edges, whenever the summand $`M_{i_r}`$ contains $`a_{rs}`$ copies of the summand $`M_{i_s^{}}`$ under the embedding $`\phi _i`$. As $`i`$ varies, one obtains an infinite graph called a Bratteli diagram of the $`AF`$-algebra . The Bratteli diagram defines a unique $`AF`$-algebra.

D. The stationary $`AF`$-algebras. If the homomorphisms $`\phi _1=\phi _2=\mathrm{}=Const`$ in the definition of the $`AF`$-algebra $`A`$, the $`AF`$-algebra $`A`$ is called stationary. The Bratteli diagram of a stationary $`AF`$-algebra looks like a periodic graph with the incidence matrix $`A=(a_{rs})`$ repeated over and over again. Since matrix $`A`$ is a non-negative integer matrix, one can take a power of $`A`$ to obtain a strictly positive integer matrix – which we always assume to be the case. The stationary $`AF`$-algebra has a non-trivial group of the automorphisms , Ch.6.

### 2.2 Teichmueller functor

A. The Jacobi-Perron continued fractions. Let $`\lambda =(\lambda _1,\mathrm{},\lambda _n)`$ be a vector with the non-negative real entries, such that $`\lambda _10`$; consider a projective class $`(1,\theta _1,\mathrm{},\theta _{n1})`$ of $`\lambda `$, where $`\theta _{i1}=\frac{\lambda _i}{\lambda _1}`$ for $`1in`$. The continued fraction

$$\left(\begin{array}{c}1\\ \theta _1\\ \mathrm{}\\ \theta _{n1}\end{array}\right)=\underset{k\mathrm{}}{lim}\left(\begin{array}{ccccc}0& 0& \mathrm{}& 0& 1\\ 1& 0& \mathrm{}& 0& b_1^{(1)}\\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 1& b_{n1}^{(1)}\end{array}\right)\mathrm{}\left(\begin{array}{ccccc}0& 0& \mathrm{}& 0& 1\\ 1& 0& \mathrm{}& 0& b_1^{(k)}\\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 1& b_{n1}^{(k)}\end{array}\right)\left(\begin{array}{c}0\\ 0\\ \mathrm{}\\ 1\end{array}\right),$$

where $`b_i^{(j)}N\{0\}`$, is called the Jacobi-Perron fraction. To recover the integers $`b_i^{(k)}`$ from the vector $`(\theta _1,\mathrm{},\theta _{n1})`$, one has to repeatedly solve the following system of equations: $`\theta _1=b_1^{(1)}+\frac{1}{\theta _{n1}^{}}`$, $`\theta _2=b_2^{(1)}+\frac{\theta _1^{}}{\theta _{n1}^{}}`$, $`\theta _{n1}=b_{n1}^{(1)}+\frac{\theta _{n2}^{}}{\theta _{n1}^{}}`$, where $`(\theta _1^{},\mathrm{},\theta _{n1}^{})`$ is the next input vector. Thus, each vector $`(\theta _1,\mathrm{},\theta _{n1})`$ gives rise to a formal Jacobi-Perron continued fraction, which can be convergent or divergent. We let $`A^{(0)}=\delta _{ij}`$ (the Kronecker delta) and $`A_i^{(k+n)}=_{j=0}^{n1}b_i^{(k)}A_i^{(\nu +j)},b_0^{(k)}=1,`$ where $`i=0,\mathrm{},n1`$ and $`k=0,1,\mathrm{},\mathrm{}`$. The Jacobi-Perron continued fraction of vector $`(\theta _1,\mathrm{},\theta _{n1})`$ is said to be convergent, if $`\theta _i=lim_k\mathrm{}\frac{A_i^{(k)}}{A_0^{(k)}}`$ for all $`i=1,\mathrm{},n1`$. Unless $`n=2`$, the convergence of the individual Jacobi-Perron fraction is a difficult problem; however, it is known that the Jacobi-Perron fractions converge for a generic subset of the vectors in the space $`R^{n1}`$ .

B. The toric $`AF`$-algebras. Denote by $`T_S(g)`$ the Teichmueller space of genus $`g1`$ with a distinguished point $`S`$. Let $`qH^0(S,\mathrm{\Omega }^2)`$ be a holomorphic quadratic differential on the Riemann surface $`S`$, such that all zeroes of $`q`$ (if any) are simple. By $`\stackrel{~}{S}`$ we mean a double cover of $`S`$ ramified over the zeroes of $`q`$ and by $`H_1^{odd}(\stackrel{~}{S})`$ the odd part of the integral homology of $`\stackrel{~}{S}`$ relative to the zeroes. Note that $`H_1^{odd}(\stackrel{~}{S})Z^n`$, where $`n=6g6`$ if $`g2`$ and $`n=2`$ if $`g=1`$. It follows from the Main Theorem of , that $`T_S(g)\{pt\}Hom(H_1^{odd}(\stackrel{~}{S});R)\{0\}`$, where $`0`$ is the zero homomorphism 222To be precise, the theorem mentions a local homeomorphism $`h:Hom(H_1^{odd}(\stackrel{~}{S});R)\{0\}T_S(g)\{pt\}`$ , p.222. Since $`T_S(g)`$ is simply connected, $`h`$ extends to a global homeomorphism between the two spaces. Indeed, let $`\lambda Hom(H_1^{odd}(\stackrel{~}{S});R)`$. It is easy to see, that $`\lambda =0`$ corresponds to the distinguished point $`ST_S(g)`$, while $`\lambda =\mathrm{}`$ represent the boundary of the space $`T_S(g)`$; thus, every ball $`|\lambda |<C`$ is homotopy equivalent to the ball $`|\lambda |<\mathrm{}`$. Note also, that we are interested in $`q`$’s with the generic (simple) zeroes; the higher order zeroes – which can be an obstacle in the construction of global coordinates – are excluded. The interested reader can consult , p.425 (the last paragraph) for the details; the mentioned there piecewise linear integral structure breaks $`T_S(g)`$ into a finite number of cones issued from $`S`$ and in these terms our construction means that we take a cone and extend it (by the linearity) to the entire $`T_S(g)`$. The cones differ from each other by a permutation on the set $`\lambda =(\lambda _1,\mathrm{},\lambda _n)`$; distinct permutations correspond to the different coordinates in the space $`T_S(g)`$. . Finally, denote by $`\lambda =(\lambda _1,\mathrm{},\lambda _n)`$ the image of a basis of $`H_1^{odd}(\stackrel{~}{S})`$ in the real line $`R`$, such that $`\lambda _10`$. (Note that such an option always exists, since the zero homomorphism is excluded.) We let $`\theta =(\theta _1,\mathrm{},\theta _{n1})`$, where $`\theta _i=\lambda _{i1}/\lambda _1`$. Recall that, up to a scalar multiple, the vector $`(1,\theta )R^n`$ is the limit of a generically convergent Jacobi-Perron continued fraction:

$$\left(\begin{array}{c}1\\ \theta \end{array}\right)=\underset{k\mathrm{}}{lim}\left(\begin{array}{cc}0& 1\\ I& b_1\end{array}\right)\mathrm{}\left(\begin{array}{cc}0& 1\\ I& b_k\end{array}\right)\left(\begin{array}{c}0\\ I\end{array}\right),$$

where $`b_i=(b_1^{(i)},\mathrm{},b_{n1}^{(i)})^T`$ is a vector of the non-negative integers, $`I`$ the unit matrix and $`I=(0,\mathrm{},0,1)^T`$. We introduce an $`AF`$-algebra, $`A_\theta `$, via the Bratteli diagram, shown in Fig.1. (The numbers $`b_j^{(i)}`$ of the diagram indicate the multiplicity of edges of the graph.) Let us call $`A_\theta `$ a toric $`AF`$-algebra. Note that in the $`g=1`$ case, the Jacobi-Perron fraction coincides with the regular continued fraction and $`A_\theta `$ becomes the Effros-Shen $`AF`$-algebra of a noncommutative torus .

C. The Teichmueller functor. Denote by $`V`$ the maximal subset of $`T_S(g)`$, such that for each Riemann surface $`RV`$, there exists a convergent Jacobi-Perron continued fraction. Let $`F`$ be the map which sends the Riemann surfaces into the toric $`AF`$-algebras according to the formula $`RA_\theta `$. Finally, let $`W`$ be the image of $`V`$ under the mapping $`F`$.

###### Lemma 1

() The set $`V`$ is a generic subset of $`T_S(g)`$ and the map $`F`$ has the following properties:

(i) $`VW\times (0,\mathrm{})`$ is a trivial fiber bundle, whose projection map $`p:VW`$ coincides with $`F`$;

(ii) $`F:VW`$ is a covariant functor, which maps isomorphic Riemann surfaces $`R,R^{}V`$ to stably isomorphic toric $`AF`$-algebras $`A_\theta ,A_\theta ^{}W`$.

## 3 Proof of theorem 1

As before, let $`W`$ denote the set of toric $`AF`$-algebras of genus $`g2`$. Let $`G`$ be a finitely presented group and $`G\times WW`$ its action on $`W`$ by the stable isomorphisms of toric $`AF`$-algebras; in other words, $`\gamma (A_\theta )KA_\theta K`$ for all $`\gamma G`$ and all $`A_\theta W`$. The following preparatory lemma will be important.

###### Lemma 2

For each $`A_\theta W`$, there exists a representation $`\rho _{A_\theta }:GGL_{6g6}(Z)`$.

Proof. The proof of lemma is based on the following well known criterion of the stable isomorphism for the (toric) $`AF`$-algebras: a pair of such algebras $`A_\theta ,A_\theta ^{}`$ are stably isomorphic if and only if their Bratteli diagrams coincide, except (possibly) a finite part of the diagram, see , Theorem 2.3. (Note, that the order isomorphism between the dimension groups, mentioned in the original text, translates to the language of the Bratteli diagrams as stated.)

Let $`G`$ be a finitely presented group on the generators $`\{\gamma _1,\mathrm{},\gamma _m\}`$ subject to relations $`r_1,\mathrm{},r_n`$. Let $`A_\theta W`$. Since $`G`$ acts on the toric $`AF`$-algebra $`A_\theta `$ by stable isomorphisms, the toric $`AF`$-algebras $`A_{\theta _1}:=\gamma _1(A_\theta ),\mathrm{},A_{\theta _m}:=\gamma _m(A_\theta )`$ are stably isomorphic to $`A_\theta `$; moreover, by transitivity, they are also pairwise stably isomorphic. Therefore, the Bratteli diagrams of $`A_{\theta _1},\mathrm{},A_{\theta _m}`$ coincide everywhere except, possibly, some finite parts. We shall denote by $`A_{\theta _{\mathrm{max}}}W`$ a toric $`AF`$-algebra, whose Bratteli diagram is the maximal common part of the Bratteli diagrams of $`A_{\theta _i}`$ for $`1im`$; such a choice is unique and defined correctly because the set $`\{A_{\theta _i}\}`$ is a finite set. By the definition of a toric $`AF`$-algebra, the vectors $`\theta _i=(1,\theta _1^{(i)},\mathrm{},\theta _{6g7}^{(i)})`$ are related to the vector $`\theta _{\mathrm{max}}=(1,\theta _1^{(\mathrm{max})},\mathrm{},\theta _{6g7}^{(\mathrm{max})})`$ by the formula:

$$\left(\begin{array}{c}1\\ \theta _1^{(i)}\\ \mathrm{}\\ \theta _{6g7}^{(i)}\end{array}\right)=\underset{A_i}{\underset{}{\left(\begin{array}{ccccc}0& 0& \mathrm{}& 0& 1\\ 1& 0& \mathrm{}& 0& b_1^{(1)(i)}\\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 1& b_{6g7}^{(1)(i)}\end{array}\right)\mathrm{}\left(\begin{array}{ccccc}0& 0& \mathrm{}& 0& 1\\ 1& 0& \mathrm{}& 0& b_1^{(k)(i)}\\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 1& b_{6g7}^{(k)(i)}\end{array}\right)}}\left(\begin{array}{c}1\\ \theta _1^{(\mathrm{max})}\\ \mathrm{}\\ \theta _{6g7}^{(\mathrm{max})}\end{array}\right)$$

The above expression can be written in the matrix form $`\theta _i=A_i\theta _{\mathrm{max}}`$, where $`A_iGL_{6g6}(Z)`$. Thus, one gets a matrix representation of the generator $`\gamma _i`$, given by the formula $`\rho _{A_\theta }(\gamma _i):=A_i`$.

The map $`\rho _{A_\theta }:GGL_{6g6}(Z)`$ extends to the rest of the group $`G`$ via its values on the generators; namely, for every $`gG`$ one sets $`\rho _{A_\theta }(g)=A_1^{k_1}\mathrm{}A_m^{k_m}`$, whenever $`g=\gamma _1^{k_1}\mathrm{}\gamma _m^{k_m}`$. Let us verify, that the map $`\rho _{A_\theta }`$ is a well defined homomorphism of groups $`G`$ and $`GL_{6g6}(Z)`$. Indeed, let us write $`g_1=\gamma _1^{k_1}\mathrm{}\gamma _m^{k_m}`$ and $`g_2=\gamma _1^{s_1}\mathrm{}\gamma _m^{s_m}`$ for a pair of elements $`g_1,g_2G`$; then their product $`g_1g_2=\gamma _1^{k_1}\mathrm{}\gamma _m^{k_m}\gamma _1^{s_1}\mathrm{}\gamma _m^{s_m}=\gamma _1^{l_1}\mathrm{}\gamma _m^{l_m}`$, where the last equality is obtained by a reduction of words using the relations $`r_1,\mathrm{},r_n`$. One can write relations $`r_i`$ in their matrix form $`\rho _{A_\theta }(r_i)`$; thus, one gets the matrix equality $`A_1^{l_1}\mathrm{}A_m^{l_m}=A_1^{k_1}\mathrm{}A_m^{k_m}A_1^{s_1}\mathrm{}A_m^{s_m}`$. It is immediate from the last equation, that $`\rho _{A_\theta }(g_1g_2)=A_1^{l_1}\mathrm{}A_m^{l_m}=A_1^{k_1}\mathrm{}A_m^{k_m}A_1^{s_1}\mathrm{}A_m^{s_m}=\rho _{A_\theta }(g_1)\rho _{A_\theta }(g_2)`$ for $`g_1,g_2G`$, i.e. $`\rho _{A_\theta }`$ is a homomorphism. Lemma 2 follows. $`\mathrm{}`$

Let $`W_{aper}W`$ be a set consisting of the toric $`AF`$-algebras, whose Bratteli diagrams do not contain periodic (infinitely repeated) blocks; these are known as non-stationary toric $`AF`$-algebras and they are generic in the set $`W`$ endowed with the natural topology. We call the action of $`G`$ on the toric $`AF`$-algebra $`A_\theta W`$ free, if $`\gamma (A_\theta )=A_\theta `$ implies $`\gamma =Id`$.

###### Lemma 3

Let $`A_\theta W_{aper}`$ and $`G`$ be free on the $`A_\theta `$. Then $`\rho _{A_\theta }`$ is a faithful representation.

Proof. Since the action of $`G`$ is free, to prove that $`\rho _{A_\theta }`$ is faithful, it remains to show, that in the formula $`\theta _i=A_i\theta _{\mathrm{max}}`$, it holds $`A_i=I`$, if and only if, $`\theta _i=\theta _{\mathrm{max}}`$, where $`I`$ is the unit matrix. Indeed, it is immediate that $`A_i=I`$ implies $`\theta _i=\theta _{\mathrm{max}}`$. Suppose now that $`\theta _i=\theta _{\mathrm{max}}`$ and, let to the contrary, $`A_iI`$. One gets $`\theta _i=A_i\theta _{\mathrm{max}}=\theta _{\mathrm{max}}`$. Such an equation has a non-trivial solution, if and only if, the vector $`\theta _{\mathrm{max}}`$ has a periodic Jacobi-Perron fraction; the period of such a fraction is given by the matrix $`A_i`$. This is impossible, since it has been assumed, that $`A_{\theta _{\mathrm{max}}}W_{aper}`$. The contradiction finishes the proof of lemma 3. $`\mathrm{}`$

Let $`G=Mod(X)`$, where $`X`$ is a surface of genus $`g2`$. The group $`G`$ is finitely presented ; it acts on the Teichmueller space $`T(g)`$ by isomorphisms of the Riemann surfaces. Moreover, the action of $`G`$ is free on a generic set, $`UT(g)`$, consisting of the Riemann surfaces with the trivial group of the automorphisms. Recall, that there exists a functor $`F:VW`$ between the Riemann surfaces and toric $`AF`$-algebras, see lemma 1.

###### Lemma 4

The pre-image $`F^1(W_{aper})`$ is a generic set in the space $`T(g)`$.

Proof. Note, that the set of stationary toric $`AF`$-algebras is a countable set. The functor $`F`$ is a surjective map, which is continuous with respect to the natural topology on the sets $`V`$ and $`W`$. Therefore, the pre-image of the complement of a countable set is a generic set. $`\mathrm{}`$

To finish the proof, consider the set $`UF^1(W_{aper})`$; this is a non-empty set, since it is the intersection of two generic subsets of $`T(g)`$. Let $`R`$ be a point (a Riemann surface) in the above set. In view of lemma 1, group $`G`$ acts on the toric $`AF`$-algebra $`A_\theta =F(R)`$ by the stable isomorphisms. By the construction, the action is free and $`A_\theta W_{aper}`$. In view of lemma 3, one gets a faithful representation $`\rho =\rho _{A_\theta }`$ of the group $`G=Mod(X)`$ into the matrix group $`GL_{6g6}(Z)`$. Theorem 1 follows. $`\mathrm{}`$

The Fields Institute for Mathematical Sciences, Toronto, ON, Canada, E-mail: igor.v.nikolaev@gmail.com

Current address: 101-315 Holmwood Ave., Ottawa, ON, Canada, K1S 2R2


format: text

text:
# 1 Introduction

## 1 Introduction

In a landmark paper in 1991, Jackson and Fersht  demonstrated that chymotrypsin inhibitor 2 folds without significantly populating any meta-stable intermediate state. Since then, it has become clear that this protein is far from unique; the same behavior has been observed for many small single-domain proteins . It is tempting to interpret the apparent two-state behavior of these proteins in terms of a simple free-energy landscape with two minima separated by a single barrier, where the minima represent the native and unfolded states, respectively. If the barrier is high, this picture provides an explanation of why the folding kinetics are single exponential, and why the folding thermodynamics show two-state character.

However, it is well-known that the free-energy barrier, $`\mathrm{\Delta }F`$, is not high for all these proteins. In fact, assuming the folding time $`\tau _\text{f}`$ to be given by $`\tau _\text{f}=\tau _0\mathrm{exp}(\mathrm{\Delta }F/kT)`$ with $`\tau _01\mu s`$ , it is easy to find examples of proteins with $`\mathrm{\Delta }F`$ values of a few $`kT`$  ($`k`$ is Boltzmann’s constant and $`T`$ the temperature). It should also be mentioned that Garcia-Mira et al.  recently found a protein that appears to fold without crossing any free-energy barrier.

Suppose the native and unfolded states coexist at the folding temperature and that there is no well-defined intermediate state, but that a clear free-energy barrier is missing. What type of folding behavior should one then expect? In particular, would such a protein, due to the lack of a clear free-energy barrier, show easily detectable deviations from two-state thermodynamics and single-exponential kinetics? Here we investigate this question based on Monte Carlo simulations of a designed three-helix-bundle protein .

Our study consists of three parts. First, we investigate whether or not melting curves for this model protein show two-state character. Second, we ask whether the relaxation behavior is single exponential or not, based on ensemble kinetics at the folding temperature. Third, inspired by energy-landscape theory (for a recent review, see Refs. ), we try to interpret the folding dynamics of this system in terms of simple diffusive motion in a low-dimensional free-energy landscape.

## 2 Model and Methods

### 2.1 The Model

Simulating atomic models for protein folding remains a challenge, although progress is currently being made in this area . Here, for computational efficiency, we consider a reduced model with 5–6 atoms per amino acid , in which the side chains are replaced by large Cβ atoms. Using this model, we study a designed three-helix-bundle protein with 54 amino acids.

The model has the Ramachandran torsion angles $`\varphi _i`$, $`\psi _i`$ as its degrees of freedom, and is sequence-based with three amino acid types: hydrophobic (H), polar (P) and glycine (G). The sequence studied consists of three identical H/P segments with 16 amino acids each (PPHPPHHPPHPPHHPP), separated by two short GGG segments . The H/P segment is such that it can make an $`\alpha `$-helix with all the hydrophobic amino acids on the same side.

The interaction potential

$$E=E_{\text{loc}}+E_{\text{ev}}+E_{\text{hb}}+E_{\text{hp}}$$
(1)
is composed of four terms. The local potential $`E_{\text{loc}}`$ has a standard form with threefold symmetry,

$$E_{\text{loc}}=\frac{ϵ_\varphi }{2}\underset{i}{}(1+\mathrm{cos}3\varphi _i)+\frac{ϵ_\psi }{2}\underset{i}{}(1+\mathrm{cos}3\psi _i).$$
(2)
The excluded-volume term $`E_{\text{ev}}`$ is given by a hard-sphere potential of the form

$$E_{\text{ev}}=ϵ_{\text{ev}}\stackrel{}{}_{i<j}\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12},$$
(3)
where the sum runs over all possible atom pairs except those consisting of two hydrophobic Cβ. The parameter $`\sigma _{ij}`$ is given by $`\sigma _{ij}=\sigma _i+\sigma _j+\mathrm{\Delta }\sigma _{ij}`$, where $`\mathrm{\Delta }\sigma _{ij}=0.625`$ Å for CβC, CβN and CβO pairs that are connected by a sequence of three covalent bonds, and $`\mathrm{\Delta }\sigma _{ij}=0`$ Å otherwise. The introduction of the parameter $`\mathrm{\Delta }\sigma _{ij}`$ can be thought of as a change of the local potential.

The hydrogen-bond term $`E_{\text{hb}}`$ has the form

$$E_{\text{hb}}=ϵ_{\text{hb}}\underset{ij}{}u(r_{ij})v(\alpha _{ij},\beta _{ij}),$$
(4)
where the functions $`u(r)`$ and $`v(\alpha ,\beta )`$ are given by

$`u(r)`$ $`=`$ $`5\left({\displaystyle \frac{\sigma _{\text{hb}}}{r}}\right)^{12}6\left({\displaystyle \frac{\sigma _{\text{hb}}}{r}}\right)^{10}`$ (5)
$`v(\alpha ,\beta )`$ $`=`$ $`\{\begin{array}{cc}\mathrm{cos}^2\alpha \mathrm{cos}^2\beta \hfill & \alpha ,\beta >90^{}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (8)
The sum in Eq. 4 runs over all possible HO pairs, and $`r_{ij}`$ denotes the HO distance, $`\alpha _{ij}`$ the NHO angle, and $`\beta _{ij}`$ the HOC angle. The last term of the potential, the hydrophobicity term $`E_{\text{hp}}`$, is given by

$$E_{\text{hp}}=ϵ_{\text{hp}}\underset{i<j}{}\left[\left(\frac{\sigma _{\text{hp}}}{r_{ij}}\right)^{12}2\left(\frac{\sigma _{\text{hp}}}{r_{ij}}\right)^6\right],$$
(9)
where the sum runs over all pairs of hydrophobic Cβ.

To speed up the calculations, a cutoff radius $`r_c`$ is used, which is taken to be 4.5 Å for $`E_{\text{ev}}`$ and $`E_{\text{hb}}`$, and 8 Å for $`E_{\text{hp}}`$. Numerical values of all energy and geometry parameters can be found elsewhere .

The thermodynamic behavior of this three-helix-bundle protein has been studied before . These studies demonstrated that this model protein has the following properties:

* It does form a stable three-helix bundle, except for a twofold topological degeneracy. These two topologically distinct states both contain three right-handed helices. They differ in how the helices are arranged. If we let the first two helices form a U, then the third helix is in front of the U in one case (FU), and behind the U in the other case (BU). The reason that the model is unable to discriminate between these two states is that their contact maps are effectively very similar .
* It makes more stable helices than the corresponding one- and two-helix sequences, which is in accord with the experimental fact that tertiary interactions generally are needed for secondary structure to become stable.
* It undergoes a first-order-like folding transition directly from an expanded state to the three-helix-bundle state, without any detectable intermediate state. At the folding temperature $`T_\text{f}`$, there is a pronounced peak in the specific heat.

Here we analyze the folding dynamics of this protein in more detail, through an extended study of both thermodynamics and kinetics.

As a measure of structural similarity with the native state, we monitor a parameter $`Q`$ that we call nativeness (the same as in ). To calculate $`Q`$, we use representative conformations for the FU and BU topologies, respectively, obtained by energy minimization. For a given conformation, we compute the root-mean-square deviations $`\delta _{\text{FU}}`$ and $`\delta _{\text{BU}}`$ from these two representative conformations (calculated over all backbone atoms). The nativeness $`Q`$ is then obtained as

$$Q=\mathrm{max}[\mathrm{exp}\left(\delta _{\text{FU}}^2/(10\text{Å})^2\right),\mathrm{exp}\left(\delta _{\text{BU}}^2/(10\text{Å})^2\right)],$$
(10)
which makes $`Q`$ a dimensionless number between 0 and 1.

Energies are quoted in units of $`kT_\text{f}`$, with the folding temperature $`T_\text{f}`$ defined as the specific heat maximum. In the dimensionless energy unit used in our previous study , this temperature is given by $`kT_\text{f}=0.6585\pm 0.0006`$.

### 2.2 Monte Carlo Methods

To simulate the thermodynamic behavior of this model, we use simulated tempering , in which the temperature is a dynamic variable. This method is chosen in order to speed up the calculations at low temperatures. Our simulations are started from random configurations. The temperatures studied range from $`0.95T_\text{f}`$ to $`1.37T_\text{f}`$.

The temperature update is a standard Metropolis step. In conformation space we use two different elementary moves: first, the pivot move in which a single torsion angle is turned; and second, a semi-local method  that works with seven or eight adjacent torsion angles, which are turned in a coordinated manner. The non-local pivot move is included in our calculations in order to accelerate the evolution of the system at high temperatures.

Our kinetic simulations are also Monte Carlo-based, and only meant to mimic the time evolution of the system in a qualitative sense. They differ from our thermodynamic simulations in two ways: first, the temperature is held constant; and second, the non-local pivot update is not used, but only the semi-local method . This restriction is needed in order to avoid large unphysical deformations of the chain.

Statistical errors on thermodynamic results are obtained by jackknife analysis  of results from ten or more independent runs, each containing several folding/unfolding events. All errors quoted are 1$`\sigma `$ errors. Statistical errors on relaxation times are difficult to determine due to uncertainties about where the large-time behavior sets in and are therefore omitted. We estimate that the uncertainties on our calculated relaxation times are about 10%. The statistical errors on the results obtained by numerical solution of the diffusion equation are, however, significantly smaller than this.

All fits of data discussed below are carried out by using a Levenberg-Marquardt procedure .

### 2.3 Analysis

Melting curves for proteins are often described in terms of a two-state picture. In the two-state approximation, the average of a quantity $`X`$ at temperature $`T`$ is given by

$$X(T)=\frac{X_\text{u}+X_\text{n}K(T)}{1+K(T)},$$
(11)
where $`K(T)=P_\text{n}(T)/P_\text{u}(T)`$, $`P_\text{n}(T)`$ and $`P_\text{u}(T)`$ being the populations of the native and unfolded states, respectively. Likewise, $`X_\text{n}`$ and $`X_\text{u}`$ denote the respective values of $`X`$ in the native and unfolded states. The effective equilibrium constant $`K(T)`$ is to leading order given by $`K(T)=\mathrm{exp}[(1/kT1/kT_\text{m})\mathrm{\Delta }E]`$, where $`T_\text{m}`$ is the midpoint temperature and $`\mathrm{\Delta }E`$ the energy difference between the two states. With this $`K(T)`$, a fit to Eq. 11 has four parameters: $`\mathrm{\Delta }E`$, $`T_\text{m}`$ and the two baselines $`X_\text{u}`$ and $`X_\text{n}`$.

A simple but powerful method for quantitative analysis of the folding dynamics is obtained by assuming the motion along different reaction coordinates to be diffusive . The folding process is then modeled as one-dimensional Brownian motion in an external potential given by the free energy $`F(r)=kT\mathrm{ln}P_{\text{eq}}(r)`$, where $`P_{\text{eq}}(r)`$ denotes the equilibrium distribution of $`r`$. Thus, it is assumed that the probability distribution of $`r`$ at time $`t`$, $`P(r,t)`$, obeys Smoluchowski’s diffusion equation

$$\frac{P(r,t)}{t}=\frac{}{r}\left[D(r)\left(\frac{P(r,t)}{r}+\frac{P(r,t)}{kT}\frac{F(r)}{r}\right)\right],$$
(12)
where $`D(r)`$ is the diffusion coefficient.

This picture is not expected to hold on short timescales, due to the projection onto a single coordinate $`r`$, but may still be useful provided that the diffusive behavior sets in on a timescale that is small compared to the relaxation time. By estimating $`D(r)`$ and $`F(r)`$, it is then possible to predict the relaxation time from Eq. 12. Such an analysis has been successfully carried through for a lattice protein .

The relaxation behavior predicted by Eq. 12 is well understood when $`F(r)`$ has the shape of a double well with a clear barrier. In this situation, the relaxation is single exponential with a rate constant given by Kramers’ well-known result . However, this result cannot be applied to our model, in which the free-energy barrier is small or absent, depending on which reaction coordinate is used. Therefore, we perform a detailed study of Eq. 12 for some relevant choices of $`D(r)`$ and $`F(r)`$, using analytical as well as numerical methods.

## 3 Results

### 3.1 Thermodynamics

In our thermodynamic analysis, we study the five different quantities listed in Table 1. The first question we ask is to what extent the temperature dependence of these quantities can be described in terms of a first-order two-state system (see Eq. 11).

Fits of our data to this equation show that the simple two-state picture is not perfect ($`\chi ^2`$ per degree of freedom, dof, of $``$ 10), but this can be detected only because the statistical errors are very small at high temperatures ($`<0.1\%`$). In fact, if we assign artifical statistical errors of 1% to our data points, an error size that is not uncommon for experimental data, then the fits become perfect with a $`\chi ^2`$/dof close to unity. Fig. 1 shows the temperature dependence of the hydrogen-bond energy $`E_{\text{hb}}`$ and the radius of gyration $`R_\text{g}`$, along with our two-state fits.

Table 1 gives a summary of our two-state fits. In particular, we see that the fitted values of both the energy change $`\mathrm{\Delta }E`$ and the midpoint temperature $`T_\text{m}`$ are similar for the different quantities. It is also worth noting that the $`T_\text{m}`$ values fall close to the folding temperature $`T_\text{f}`$, defined as the maximum of the specific heat. The difference between the highest and lowest values of $`T_\text{m}`$ is less than 1%. There is a somewhat larger spread in $`\mathrm{\Delta }E`$, but this parameter has a larger statistical error.

So, the melting curves show two-state character, and the fitted parameters $`\mathrm{\Delta }E`$ and $`T_\text{m}`$ are similar for different quantities. From this it may be tempting to conclude that the thermodynamic behavior of this protein can be fully understood in terms of a two-state system. The two-state picture is, nevertheless, an oversimplification, as can be seen from the shapes of the free-energy profiles $`F(E)`$ and $`F(Q)`$. Fig. 2 shows these profiles at $`T=T_\text{f}`$. First of all, these profiles show that the native and unfolded states coexist at $`T=T_\text{f}`$, so the folding transition is first-order-like. However, there is no clear free-energy barrier separating the two states; $`F(Q)`$ exhibits a very weak barrier, $`<1kT`$, whereas $`F(E)`$ shows no barrier at all. In fact, $`F(E)`$ has the shape of a square well rather than a double well.

Phase transition terminology is, by necessity, ambiguous for a finite system like this, but if states with markedly different $`E`$ or $`Q`$ coexist it does make sense to call the transition first-order-like, even if a free-energy barrier is missing. At a second-order phase transition, the free-energy profile is wide, but the minimum remains unique.

### 3.2 Kinetics

Our kinetic study is performed at $`T=T_\text{f}`$. Using Monte Carlo dynamics (see Model and Methods), we study the relaxation of ensemble averages of various quantities. For this purpose, we performed a set of 3000 folding simulations, starting from equilibrium conformations at temperature $`T_01.06T_\text{f}`$. At this temperature, the chain is extended and has a relatively low secondary-structure content (see Fig. 1).

In the absence of a clear free-energy barrier (see Fig. 2), it is not obvious whether or not the relaxation should be single exponential. To get an idea of what to expect for a system like this, we consider the relaxation of the energy $`E`$ in a potential $`F(E)`$ that has the form of a perfect square well at $`T=T_\text{f}`$. For this idealized $`F(E)`$ and a constant diffusion coefficient $`D(E)`$, it is possible to solve Eq. 12 analytically for relaxation at an arbitrary temperature $`T`$. This solution is given in Appendix A, for the initial condition that $`P(E,t=0)`$ is the equilibrium distribution at temperature $`T_0`$. Using this result, the deviation from single-exponential behavior can be mapped out as a function of $`T_0`$ and $`T`$, as is illustrated in Fig. 3. The size of the deviation depends on both $`T_0`$ and $`T`$, but is found to be small for a wide range of $`T_0`$, $`T`$ values. This clearly demonstrates that the existence of a free-energy barrier is not a prerequisite to observe single-exponential relaxation.

Let us now turn to the results of our simulations. Fig. 4 shows the relaxation of the average energy $`E`$ and the average nativeness $`Q`$ in Monte Carlo (MC) time. In both cases, the large-time data can be fitted to a single exponential, which gives relaxation times of $`\tau 1.710^7`$ and $`\tau 1.810^7`$ for $`E`$ and $`Q`$, respectively, in units of elementary MC steps. The corresponding fits for the radius of gyration and the hydrogen-bond energy (data not shown) give relaxation times of $`\tau 2.110^7`$ and $`\tau 1.810^7`$, respectively. The fit for the radius of gyration has a larger uncertainty than the others, because the data points have larger errors in this case.

The differences between our four fitted $`\tau `$ values are small and most probably due to limited statistics for the large-time behavior. Averaging over the four different variables, we obtain a relaxation time of $`\tau 1.810^7`$ MC steps for this protein. The fact that the relaxation times for the hydrogen-bond energy and the radius of gyration are approximately the same shows that helix formation and chain collapse proceed in parallel for this protein. This finding is in nice agreement with recent experimental results for small helical proteins .

For $`Q`$, it is necessary to go to very short times in order to see any significant deviation from a single exponential (see Fig. 4). For $`E`$, we find that the single-exponential behavior sets in at roughly $`\tau /3`$, which means that the deviation from this behavior is larger than in the analytical calculation above. On the other hand, for comparisons with experimental data, we expect the behavior of $`Q`$ to be more relevant than that of $`E`$. The simulations confirm that the relaxation can be approximately single exponential even if there is no clear free-energy barrier.

To translate the relaxation time for this protein into physical units, we compare with the reconfiguration time for the corresponding one-helix segment. To that end, we performed a kinetic simulation of this 16-amino acid segment at the same temperature, $`T=T_\text{f}`$. This temperature is above the midpoint temperature for the one-helix segment, which is $`0.95T_\text{f}`$ . So, the isolated one-helix segment is unstable at $`T=T_\text{f}`$, but makes frequent visits to helical states with low hydrogen-bond energy, $`E_{\text{hb}}`$. To obtain the reconfiguration time, we fitted the large-time behavior of the autocorrelation function for $`E_{\text{hb}}`$,

$$C_{\text{hb}}(t)=E_{\text{hb}}(t)E_{\text{hb}}(0)E_{\text{hb}}(0)^2,$$
(13)
to an exponential. The exponential autocorrelation time, which can be viewed as a reconfiguration time, turned out to be $`\tau _\text{h}1.010^6`$ MC steps. This is roughly a factor 20 shorter than the relaxation time $`\tau `$ for the full three-helix bundle. Assuming the reconfiguration time for an individual helix to be $``$$`0.2\mu `$, we obtain relaxation and folding times of $``$$`4\mu `$s and $``$$`8\mu `$s, respectively, for the three-helix bundle. This is fast but not inconceivable for a small helical protein . In fact, the B domain of staphylococcal protein A is a three-helix-bundle protein that has been found to fold in $`<10\mu `$s, at 37.

### 3.3 Relaxation-Time Predictions

We now turn to the question of whether the observed relaxation time can be predicted based on the diffusion equation, Eq. 12. For that purpose, we need to know not only the free energy $`F(r)`$, but also the diffusion coefficient $`D(r)`$. Socci et al.  successfully performed this analysis for a lattice protein that exhibited a relatively clear free-energy barrier. Their estimate of $`D(r)`$ involved an autocorrelation time for the unfolded state. The absence of a clear barrier separating the native and unfolded states makes it necessary to take a different approach in our case.

The one-dimensional diffusion picture is not expected to hold on short timescales, but only after coarse-graining in time. A computationally convenient way to implement this coarse-graining in time is to study block averages $`b(t)`$ defined by

$$b(t)=\frac{1}{\tau _\text{b}}\underset{ts<t+\tau _\text{b}}{}r(s)t=0,\tau _\text{b},2\tau _\text{b},\mathrm{}$$
(14)
where $`\tau _\text{b}`$ is the block size and $`r`$ is the reaction coordinate considered. The diffusion coefficient can then be estimated using $`D_\text{b}(r)=(\delta b)^2/2\tau _\text{b}`$, where the numerator is the mean-square difference between two consecutive block averages, given that the first of them has the value $`r`$.

In our calculations, we use a block size of $`\tau _\text{b}=10^6`$ MC step, corresponding to the reconfiguration time $`\tau _\text{h}`$ for an individual helix. We do not expect the dynamics to be diffusive on timescales shorter than this, due to steric traps that can occur in the formation of a helix. In order for the dynamics to be diffusive, the timescale should be such that the system can escape from these traps.

Using this block size, we first make rough estimates of the relaxation times for $`E`$ and $`Q`$ based on the result in Appendix A for a square-well potential and a constant diffusion coefficient. These estimates are given by $`\tau _{\text{pred,0}}=\mathrm{\Delta }r_{\text{sw}}^2/D_\text{b}\pi ^2`$, where $`\mathrm{\Delta }r_{\text{sw}}`$ is the width of the potential and $`D_\text{b}`$ is the average diffusion coefficient. Eq. 17 in Appendix A can be applied to other observables than $`E`$. The predicted relaxation time $`\tau _{\text{pred,0}}`$ is given by $`\tau _1`$. Our estimates of $`\mathrm{\Delta }r_{\text{sw}}`$ and $`D_\text{b}`$ can be found in Table 2, along with the resulting predictions $`\tau _{\text{pred,0}}`$. We find that these simple predictions agree with the observed relaxation times $`\tau `$ within a factor of two.

We also did the same calculation for smaller block sizes, $`\tau _\text{b}=10^0`$, $`10^1,\mathrm{}`$, $`10^5`$ MC steps. This gave $`\tau _{\text{pred,0}}`$ values smaller or much smaller than the observed $`\tau `$, signaling non-diffusive dynamics. This confirms that for $`b(t)`$ to show diffusive dynamics, $`\tau _\text{b}`$ should not be smaller than the reconfiguration time for an individual helix.

Having seen the quite good results obtained by this simple calculation, we now turn to a more detailed analysis, illustrated in Fig. 5a. The block size is the same as before, $`\tau _\text{b}=10^6`$ MC steps, but the space dependence of the diffusion coefficient $`D_\text{b}(r)`$ is now taken into account, and the potential, $`F_\text{b}(r)`$, reflects the actual distribution of block averages. The potential $`F_\text{b}(r)`$, shown in Fig. 2, is not identical to that for the unblocked variables. At a first-order-like transition, we expect free-energy minima to become more pronounced when going to the block variables, provided that the block size $`\tau _\text{b}`$ is small compared with the relaxation time, because when forming the block variables one effectively integrates out fluctuations about the respective states. The results in Fig. 2 do show this tendency, although the effect is not very strong. Fig. 5b shows the diffusion coefficient $`D_\text{b}(E)`$, which is largest at intermediate values between the native and unfolded states. The behavior of $`D_\text{b}(Q)`$ (not shown) is the same in this respect. Hence, there is no sign of a kinetic bottleneck to folding for this protein.

Given $`D_\text{b}(r)`$ and $`F_\text{b}(r)`$, we solve Eq. 12 for $`P(r,t)`$ by using the finite-difference scheme in Appendix B. The initial distribution $`P(r,t=0)`$ is taken to be the same as in the kinetic simulations. We find that the mean of $`P(r,t)`$ shows single-exponential relaxation to a good approximation. An exponential fit of these data gives us a new prediction, $`\tau _{\text{pred}}`$, for the relaxation time.

From Table 2 it can be seen that the prediction obtained through this more elaborate calculation, $`\tau _{\text{pred}}`$, is not better than the previous one, $`\tau _{\text{pred,0}}`$, at least not in $`Q`$, despite that there exists a weak barrier in this coordinate (see Fig. 2b). This means that the barrier in $`Q`$ is too weak to be important for the relaxation rate. If the underlying diffusion picture, Eq. 12, had been perfect, $`\tau _{\text{pred}}`$ would have been equal to $`\tau `$, as obtained from the kinetic simulations. Our results show that this is not the case. At least in $`Q`$, there are significant deviations from the behavior predicted by this equation.

If more accurate relaxation time predictions are needed, there are different ways to proceed. One possible way is to simply increase the block size. However, for the calculation to be useful, the block size must remain small compared to the relaxation time. A more interesting possibility is to refine the simple diffusion picture defined by Eq. 12, in which, in particular, non-Markovian effects are ignored. Such effects may indeed affect folding times . Yet another possibility is to use a combination of a few different variables, perhaps $`E`$ and $`Q`$, instead of a single reaction coordinate . With a multidimensional representation of the folding process, non-Markovian effects could become smaller.

## 4 Summary and Discussion

We have analyzed the thermodynamics and kinetics of a designed three-helix-bundle protein, based on Monte Carlo calculations. We found that this model protein shows two-state behavior, in the sense that melting curves to a good approximation can be described by a simple two-state system and that the relaxation behavior is close to single exponential. A simple two-state picture is, nevertheless, an oversimplification, as the free-energy barrier separating the native and unfolded states is weak ($`1kT`$). The weakness of the barrier implies that a fitted two-state parameter such as $`\mathrm{\Delta }E`$ has no clear physical meaning, despite that the two-state fit looks good.

Reduced  and all-atom  models for small helical proteins have been studied by many other groups. Most of these studies relied on so-called Gō-type  potentials. It should therefore be pointed out that our model is sequence-based.

Using an extended version of this model that includes all atoms, we recently found similar results for two peptides, an $`\alpha `$-helix and a $`\beta `$-hairpin . Here the calculated melting curves could be directly compared with experimental data, and a reasonable quantitative agreement was found.

The smallness of the free-energy barrier prompted us to perform an analytical study of diffusion in a square-well potential. Here we studied the relaxation behavior at temperature $`T`$, starting from the equilibrium distribution at temperature $`T_0`$, for arbitrary $`T`$ and $`T_0`$. We found that this system shows a relaxation behavior that is close to single exponential for a wide range of $`T_0`$, $`T`$ values, despite the absence of a free-energy barrier. We also made relaxation-time predictions based on this square-well approximation. Here we took the diffusion coefficient to be constant. It was determined assuming the dynamics to be diffusive on timescales beyond the reconfiguration time for an individual helix. The predictions obtained this way were found to agree within a factor of two with observed relaxation times, as obtained from the kinetic simulations. So, this calculation, based on the two simplifying assumptions that the potential is a square well and that the diffusion coefficient is constant, gave quite good results. A more detailed calculation, in which these two additional assumptions were removed, did not give better results. This shows that the underlying diffusion picture leaves room for improvement.

Our kinetic study focused on the behavior at the folding temperature $`T_\text{f}`$, where the native and unfolded states, although not separated by a clear barrier, are very different, which makes the folding mechanism transparent. In particular, we found that helix formation and chain collapse could not be separated, which is in accord with experimental data by Krantz et al. . The difference between the native and unfolded states is much smaller at the lowest temperature we studied, $`0.95T_\text{f}`$, because the unfolded state is much more native-like here. Mayor et al.  recently reported experimental results on a three-helix-bundle protein, the engrailed homeodomain , including a characterization of its unfolded state. In particular, the unfolded state was found to have a high helix content. This study was performed at a temperature below $`0.95T_\text{f}`$. It would be very interesting to see what the unfolded state of this protein looks like near $`T_\text{f}`$. In our model, there is a significant decrease in helix content of the unfolded state as the temperature increases from $`0.95T_\text{f}`$ to $`T_\text{f}`$.

It is instructive to compare our results with those of Zhou and Karplus , who discussed two folding scenarios for helical proteins, based on a Gō-type Cα model. In their first scenario, folding is fast, without any obligatory intermediate, and helix formation occurs before chain collapse. In the second scenario, folding is slow with an obligatory intermediate on the folding pathway, and helix formation and chain collapse occur simultaneously. The behavior we find does not match any of these two scenarios. In our case, helix formation and chain collapse occur in parallel but folding is nevertheless fast and without any well-defined intermediate state.

### Acknowledgments

This work was in part supported by the Swedish Foundation for Strategic Research and the Swedish Research Council.

### Appendix A: Diffusion in a square well

Here we discuss Eq. 12 in the situation when the reaction coordinate $`r`$ is the energy $`E`$, and the potential $`F(E)`$ is a square well of width $`\mathrm{\Delta }E_{\text{sw}}`$ at $`T=T_\text{f}`$. This means that the equilibrium distribution is given by $`P_{\text{eq}}(E)\mathrm{exp}(\delta \beta E)`$ if $`E`$ is in the square well and $`P_{\text{eq}}(E)=0`$ otherwise, where $`\delta \beta =1/kT1/kT_\text{f}`$. Eq. 12 then becomes

$$\frac{P(E,t)}{t}=\frac{}{E}\left[D\left(\frac{P(E,t)}{E}+\delta \beta P(E,t)\right)\right].$$
(15)
For simplicity, the diffusion coefficient is assumed to be constant, $`D(E)=D`$. The initial distribution $`P(E,t=0)`$ is taken to be the equilibrium distribution at some temperature $`T_0`$, and we put $`\delta \beta _0=1/kT_01/kT_\text{f}`$.

By separation of variables, it is possible to solve Eq. 15 with this initial condition analytically for arbitrary values of the initial and final temperatures $`T_0`$ and $`T`$, respectively. In particular, this solution gives us the relaxation behavior of the average energy. The average energy at time $`t`$, $`E(t)`$, can be expressed in the form

$$E(t)=E+\underset{k=1}{\overset{\mathrm{}}{}}A_k\mathrm{e}^{t/\tau _k},$$
(16)
where $`E`$ denotes the equilibrium average at temperature $`T`$. A straightforward calculation shows that the decay constants in this equation are given by

$$1/\tau _k=\frac{D}{\mathrm{\Delta }E_{\text{sw}}^2}\left(\pi ^2k^2+\frac{1}{4}\delta \beta ^2\mathrm{\Delta }E_{\text{sw}}^2\right)$$
(17)
and the expansion coefficients by

$$A_k=B_k\mathrm{\Delta }E_{\text{sw}}\frac{\pi ^2k^2\left(\delta \beta \delta \beta _0\right)\mathrm{\Delta }E_{\text{sw}}}{\left(\pi ^2k^2+(\delta \beta _0\frac{1}{2}\delta \beta )^2\mathrm{\Delta }E_{\text{sw}}^2\right)\left(\pi ^2k^2+\frac{1}{4}\delta \beta ^2\mathrm{\Delta }E_{\text{sw}}^2\right)^2},$$
(18)
where

$$B_k=\frac{4\delta \beta _0\mathrm{\Delta }E_{\text{sw}}}{\mathrm{sinh}\frac{1}{2}\delta \beta _0\mathrm{\Delta }E_{\text{sw}}}\times \{\begin{array}{cc}\mathrm{cosh}\left(\frac{1}{2}(\delta \beta _0\frac{1}{2}\delta \beta )\mathrm{\Delta }E_{\text{sw}}\right)\mathrm{cosh}\frac{1}{4}\delta \beta \mathrm{\Delta }E_{\text{sw}}\hfill & \mathrm{if}k\mathrm{odd}\hfill \\ \mathrm{sinh}\left(\frac{1}{2}(\delta \beta _0\frac{1}{2}\delta \beta )\mathrm{\Delta }E_{\text{sw}}\right)\mathrm{sinh}\frac{1}{4}\delta \beta \mathrm{\Delta }E_{\text{sw}}\hfill & \mathrm{if}k\mathrm{even}\hfill \end{array}$$
(19)
Finally, the equilibrium average is

$$E=\frac{E_\text{n}+E_\text{u}}{2}+\frac{1}{\delta \beta }\frac{\mathrm{\Delta }E_{\text{sw}}}{2}\mathrm{coth}\frac{1}{2}\delta \beta \mathrm{\Delta }E_{\text{sw}},$$
(20)
where $`E_\text{n}`$ and $`E_\text{u}`$ are the lower and upper edges of the square well, respectively.

It is instructive to consider the behavior of this solution when $`|\delta \beta \delta \beta _0|1/\mathrm{\Delta }E_{\text{sw}}`$. The expression for the expansion coefficients can then be simplified to

$$A_kB_k\mathrm{\Delta }E_{\text{sw}}\frac{\pi ^2k^2(\delta \beta \delta \beta _0)\mathrm{\Delta }E_{\text{sw}}}{\left(\pi ^2k^2+\frac{1}{4}\delta \beta ^2\mathrm{\Delta }E_{\text{sw}}^2\right)^3}$$
(21)
with

$$B_k\frac{4\delta \beta _0\mathrm{\Delta }E_{\text{sw}}}{\mathrm{sinh}\frac{1}{2}\delta \beta _0\mathrm{\Delta }E_{\text{sw}}}\times \{\begin{array}{cc}\mathrm{cosh}^2\frac{1}{4}\delta \beta \mathrm{\Delta }E_{\text{sw}}\hfill & \mathrm{if}k\mathrm{odd}\hfill \\ \mathrm{sinh}^2\frac{1}{4}\delta \beta \mathrm{\Delta }E_{\text{sw}}\hfill & \mathrm{if}k\mathrm{even}\hfill \end{array}$$
(22)
Note that $`A_k`$ scales as $`k^2`$ if $`k\frac{1}{2\pi }|\delta \beta |\mathrm{\Delta }E_{\text{sw}}`$, and as $`1/k^4`$ if $`k\frac{1}{2\pi }|\delta \beta |\mathrm{\Delta }E_{\text{sw}}`$. Note also that the last factor in $`B_k`$ suppresses $`A_k`$ for even $`k`$ if $`T`$ is close to $`T_\text{f}`$. From these two facts it follows that $`|A_1|`$ is much larger than the other $`|A_k|`$ if $`T`$ is near $`T_\text{f}`$. This makes the deviation from a single exponential small.

### Appendix B: Numerical solution of the diffusion equation

To solve Eq. 12 numerically for arbitrary $`D(r)`$ and $`F(r)`$, we choose a finite-difference scheme of Crank-Nicolson type with good stability properties. To obtain this scheme we first discretize $`r`$. Put $`r_j=j\mathrm{\Delta }r`$, $`D_j=D(r_j)`$ and $`F_j=F(r_j)`$, and let $`𝐩(t)`$ be the vector with components $`p_j(t)=P(r_j,t)`$. Approximating the RHS of Eq. 12 with suitable finite differences, we obtain

$$\frac{d𝐩}{dt}=\mathrm{𝐀𝐩}(t),$$
(23)
where A is a tridiagonal matrix given by

$`(\mathrm{𝐀𝐩}(t))_j`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }r^2}}\left[D_{j+1/2}(p_{j+1}(t)p_j(t))D_{j1/2}(p_j(t)p_{j1}(t))\right]`$ (24)
$`+`$ $`{\displaystyle \frac{1}{4kT\mathrm{\Delta }r^2}}\left[D_{j+1}p_{j+1}(t)(F_{j+2}F_j)D_{j1}p_{j1}(t)(F_jF_{j2})\right]`$
Let now $`𝐩^n=𝐩(t_n)`$, where $`t_n=n\mathrm{\Delta }t`$. By applying the trapezoidal rule for integration to Eq. 23, we obtain

$$𝐩^{n+1}𝐩^n=\frac{\mathrm{\Delta }t}{2}\left(\mathrm{𝐀𝐩}^n+\mathrm{𝐀𝐩}^{n+1}\right).$$
(25)
This equation can be used to calculate how $`P(r,t)`$ evolves with time. It can be readily solved for $`𝐩^{n+1}`$ because the matrix $`𝐀`$ is tridiagonal.


text:
# Comment on ‘Observation of Bose-Einstein Condensation of Molecules’, cond-mat/0311617

Zwierlein et al. report in their recent preprint  a value for the scattering length of weakly bound molecules that are formed out of two fermionic 6Li atoms in different internal states. They observe a Thomas-Fermi density profile typical for a condensate of molecules, but there seems to be a serious discrepancy between the measured value of the measured molecular scattering length of $`a_{\text{mol}}^{\text{exp}}=8`$ nm  at a magnetic field of 770 G and the predicted one  of $`a_{\text{mol}}^{\text{theory}}=0.6`$ $`a_{\text{atom}}=120`$ nm. In this comment, we argue that this might be an indication that they have produced in fact a new many-body state.

As a first guide to explain this discrepancy, we want to point out that the Bose-Einstein condensed molecular gas is not dilute, in the sense that $`n^{1/3}R`$ is not very small ($`R`$ is the radius of the molecule and $`n`$ the density). The reported density is $`5\times 10^{13}`$cm-3 in the center of the trap and the measured radius of the molecule $`R100`$ nm, consistent with the predicted value for the atomic scattering length $`a_{\text{atom}}=200`$ nm. This gives us $`n^{1/3}R=0.368`$. The supposition that the molecular gas behaves as if the molecules were pointlike bosons is then not necessarily fulfilled and in some sense this experiment is already in the cross-over regime between BEC and BCS and not anymore in the true BEC limit. This might explain the discrepancy between the theory and the experiment, because the theory of Petrov, Salomon and Schlyapnikov.  only deals with the true BEC limit, for which it predicts $`a_{\text{mol}}^{\text{BEC}}=0.6`$ $`a_{\text{atom}}`$.

Apart from the result of Ref.  for molecules, there is a second class of papers that discusses the cross-over between BEC and BCS starting from the BCS-side  in stead. In these papers, a BCS ansatz is made for the wave function of the many body system. This BCS-ansatz seems to give meaningful results both in the weak-coupling (BCS) and strong-coupling (BEC) limit, predicting for the molecule-molecule scattering length the value of $`a_{\text{mol}}^{\text{BCS}}=2a_{\text{atom}}`$. This would result in a molecular scatting length of 400 nm in this case, even further away from what is observed experimentally.

Insight can be gained by a closer look at the way in which the scattering length was determined experimentally, namely through a measurement of the mean-field energy . To compare the experimental result with that of the two theories ($`a_{\text{mol}}^{\text{BEC}}`$ from , vs. $`a_{\text{mol}}^{\text{BCS}}=2a_{\text{atom}}`$ from ), we have to start from the mean-field energy predicted in these theories.

As stated, in order to obtain the scattering length from the experiment, the mean-field energy is measured . On the one hand, for a weakly interacting bose-gas, this mean-field energy equals

$$\mu _{\text{mol}}^{\text{BEC}}\left[n\right]=\frac{4\pi \mathrm{}^2a_{\text{mol}}n}{m_{\text{mol}}}.$$
(1)
This is the chemical potential of the homogeneous Bose-gas with density $`n`$. In the Thomas-Fermi approximation of describing the condensate, this is the only quantity that has to be known, since one assumes that

$$\mu _{\text{mol}}^{\text{BEC}}\left[n\left(𝐫\right)\right]+V_{\text{ext}}\left(𝐫\right)=\lambda ,$$
(2)
with $`V_{ext}`$ the external trapping potential an $`\lambda `$ the chemical potential of the trapped Bose gas. This allows us the calculate all the properties, such as the density profile $`n\left(𝐫\right)`$ and the total potential energy. The only assumption on which (2) is based is that $`\mu \left[n\left(𝐫\right)\right]\mathrm{}\overline{\omega },`$ where $`\overline{\omega }`$ denotes the geometric mean of the three trapping frequencies. The reported value of the mean field energy is 100 nK and $`\overline{\omega }=6.18`$ nK, so that the Thomas-Fermi expression (2) is valid.

On the other hand, we can also calculate the chemical potential starting from a BCS-ansatz. Analytical results can be obtained both in the strong-coupling (BEC) and weak-coupling (BCS) limits. In the strong-coupling BEC-limit, it is

$$\mu _{\text{mol}}^{\text{BCS}}=2\mu _{\text{atom}}=E_\text{b}+\frac{4\pi \mathrm{}^2a_{\text{atom}}n}{m_{\text{atom}}},$$
(3)
with $`E_\text{b}`$ the molecular binding energy. This is equal to (1), up to a density-independent constant, when we set $`a_{\text{mol}}=2a_{\text{atom}}`$ and $`m_{\text{mol}}=2m_{\text{atom}}`$, which agrees with the results from Refs. . Ref.  shows that in the case of the present experiments the analytical estimation deviates from the numerical result. When we do the numerical calculation for $`a_{\text{atom}}=200`$ nm and $`n=5\times 10^{13}`$cm$`^3,`$ we get $`\mu _{\text{mol}}^{\text{BCS}}=6.4`$ $`\mu `$K, which converted to a scattering length with (1) gives $`a_{\text{mol}}^{\text{BCS}}=253`$ nm, which is also not very close to the measured value.

The measured mean-field energy of a homogeneous molecular gas could be calculated in the BCS-model, using

$$E=_0^n\mu \left[n^{}\right]𝑑n^{}=\frac{2\pi \mathrm{}^2a_{\text{mol}}n^2}{m_{\text{mol}}}.$$
(4)
Keeping in mind the fact that the BCS-result is variational, we see that we can only estimate a lower bound for $`a_{\text{mol}}`$. Apparently the experiments of Zwierlein et al. show that this estimation is quite far off and this indicates that the BCS-wave function is probably far from the true state of the system.

Drawing a similar conclusion for the prediction of Petrov et al.  is not on the same footing, because it is not clear what many-body wave-function should be used in this case of a dense molecular system ($`n^{1/3}R=0.368`$). Anyway, also in this case, if we naively use the value of the molecule-molecule scattering length and ignore the dense character of the system, we would arrive at a total mean-field energy that is far higher than the measured energy and again conclude that this is an indication that the many-body state cannot be very similar to a ‘Bose-gas of molecules’.

In our view, the experiment of Ketterle and co-workers does not necessarily need to be interpreted as indicating that the true scattering amplitude between two molecules becomes small (8 nm measured vs. 120 nm predicted). Rather it can be interpreted as indicating that the energy of the achieved state is much lower than both the (variational) energy for a BCS state and the energy of a molecular condensate. In this sense, we suggest that the experiment may point to a novel many-body state, that is neither the BCS nor the molecular BEC state, and that still has to be identified.


text:
# Self-consistent stellar dynamical tori

## 1 Introduction

Constructing self-consistent collisionless equilibrium models is a key step to understand the structure and the dynamics of stellar systems. The construction process usually starts with the assignment of a phase-space distribution function (DF) obeying the Jeans theorem (and thus a solution of the collisionless Boltzmann equation that describes such “gravitational plasma”), which leads to an expression for the density as a function of the potential; models are then calculated by solving the Poisson equation, as a non-linear, second-order partial differential equation for the potential. In the presence of special symmetries (for example, in the case of spherical models) the problem can be reduced to the study of an ordinary differential equation, and the model construction is relatively straightforward. For more general symmetries the procedure is inherently difficult. In general, solutions have to be found numerically and turn out to exist only for a finite range of the parameters appearing in the adopted DF.

Here we present the basic properties of the simplest member of a new family of DFs that are able to generate axisymmetric, truncated (i.e., finite size) stellar dynamical models characterized by toroidal shapes. The complete description of our family of DFs, together with the results of N-body simulations aimed at the study of the stability of these collisionless equilibrium configurations, will be presented in a separate paper (Ciotti, Bertin, Londrillo 2004). We recall that toroidal, self-consistent models are sometimes mentioned in the literature (e.g., see Lynden-Bell 1962; the Author expressed concerns about their stability). With our models we will investigate this matter further, also considering the possible stabilizing effect of a massive dark matter halo. Astrophysical applications will address the modeling of peanut-shaped bulges (e.g., see Aronica et al. 2003, and references therein) and the problem of the central depression recently discovered (Lauer et al. 2002) in the luminosity profiles of some elliptical galaxies hosting supermassive black holes.

## 2 The models

The particular set of axisymmetric models presented here is defined by the DF

$$f(Q,J_z^2)=f_0|J_z|^{2\alpha }(QQ_0)^\beta \mathrm{\Theta }(QQ_0),$$
(1)
where $`f_0`$ is a (dimensional) physical scale, $`QJ_z^2/(2R_\mathrm{a}^2)`$, $`\mathrm{\Theta }`$ is the Heaviside step function, and $`\alpha `$ and $`\beta `$ are two dimensionless constants; the natural coordinates are cylindrical, $`(R,z,\phi )`$. In addition, $`\mathrm{\Psi }_\mathrm{T}𝐯^2/2=\mathrm{\Psi }_\mathrm{T}(v_R^2+v_z^2+v_\phi ^2)/2`$ is the binding energy per unit mass, $`J_z=Rv_\phi `$ is the axial component of the angular momentum (per unit mass), and $`Q_0`$ and $`R_\mathrm{a}`$ are a truncation energy and an anisotropy radius, respectively. The function $`\mathrm{\Psi }_\mathrm{T}=\mathrm{\Psi }_\mathrm{T}(R,z)`$ represents the total gravitational potential; for the moment, we only require that $`\mathrm{\Psi }_\mathrm{T}\mathrm{O}(r^1)`$ for $`r\mathrm{}`$, where $`r=\sqrt{R^2+z^2}`$. We allow for an external component also, i.e., $`\mathrm{\Psi }_\mathrm{T}=\mathrm{\Psi }+\mathrm{\Psi }_{\mathrm{ext}}`$, where $`\mathrm{\Psi }`$ is the potential associated with the density distribution derived from eq. (1) as given by the following eq. (4), while $`\mathrm{\Psi }_{\mathrm{ext}}`$ is taken to be a given function (equal to zero in the fully self-consistent case). Note that we can write

$$Q=\mathrm{\Psi }_\mathrm{T}\frac{v_\mathrm{m}^2}{2}\frac{v_\phi ^2}{2}\left(1+\frac{R^2}{R_\mathrm{a}^2}\right),$$
(2)
where $`v_\mathrm{m}^2v_R^2+v_z^2`$ is the square of the meridional (or poloidal) velocity component. Some dynamical properties of the systems associated with eq. (1) are the following: $`(i)`$ no net rotation is present, and $`(ii)`$ the two components $`\sigma _R^2`$ and $`\sigma _z^2`$ of the velocity dispersion tensor are equal. In addition, note that for $`\alpha =0`$ if we formally consider $`R_\mathrm{a}\mathrm{}`$ we recover the polytropic spheres, while for integer $`\alpha `$ and $`\beta `$ the distribution function is a sum of Fricke’s (1952) DFs. The functional dependence of the system structural and dynamical properties on $`\mathrm{\Psi }_\mathrm{T}`$ is obtained by integration over $`\mathrm{\Omega }(𝐱)`$, the velocity section of phase–space over which $`f0`$ at given $`𝐱`$. The coordinates $`(v,\zeta ,\xi )`$ to be used to integrate the various quantities over the rotationally symmetric (ellipsoidal) region $`\mathrm{\Omega }(𝐱)`$ are naturally defined by

$$v_R=v\mathrm{sin}\zeta \mathrm{cos}\xi ,v_z=v\mathrm{sin}\zeta \mathrm{sin}\xi ,v_\phi =\frac{v\mathrm{cos}\zeta }{\sqrt{1+R^2/R_\mathrm{a}^2}},$$
(3)
where $`0\zeta \pi `$, $`0\xi 2\pi `$. With this choice, the ellipsoidal radius $`v`$ of the velocity section $`\mathrm{\Omega }(𝐱)`$ runs in the range $`0v\sqrt{2[\mathrm{\Psi }_\mathrm{T}(𝐱)Q_0]}`$. Integration over the two angular coordinates gives

$$\rho =4\pi f_0g(R)\frac{B(\alpha +3/2,\beta +1)2^{\alpha +1/2}}{2\alpha +1}(\mathrm{\Psi }_\mathrm{T}Q_0)^{\alpha +\beta +3/2}\mathrm{\Theta }(\mathrm{\Psi }_\mathrm{T}Q_0),$$
(4)
where $`B(x,y)=_0^1t^{x1}(1t)^{y1}𝑑t`$ is the Euler Beta function,

$$g(R)=\frac{R^{2\alpha }}{(1+R^2/R_\mathrm{a}^2)^{\alpha +1/2}},$$
(5)
and for convergence it is required that $`\alpha >1/2`$, $`\beta >1`$. Moreover,

$$\rho \sigma _\mathrm{m}^2=4\pi f_0g(R)\frac{B(\alpha +5/2,\beta +1)2^{\alpha +5/2}}{(2\alpha +1)(2\alpha +3)}(\mathrm{\Psi }_\mathrm{T}Q_0)^{\alpha +\beta +5/2}\mathrm{\Theta }(\mathrm{\Psi }_\mathrm{T}Q_0)$$
(6)
and

$$\rho \sigma _\phi ^2=\frac{4\pi f_0g(R)}{1+R^2/R_\mathrm{a}^2}\frac{B(\alpha +5/2,\beta +1)2^{\alpha +3/2}}{2\alpha +3}(\mathrm{\Psi }_\mathrm{T}Q_0)^{\alpha +\beta +5/2}\mathrm{\Theta }(\mathrm{\Psi }_\mathrm{T}Q_0).$$
(7)
The number of free parameters (excluding those associated with the external potential) needed to describe a model completely is five: the three dimensional quantities $`f_0`$, $`R_\mathrm{a}`$, and $`Q_0`$, and the two dimensionless constants $`\alpha `$ and $`\beta `$. The relation between the meridional velocity dispersion $`\sigma _\mathrm{m}^2`$ and the azimuthal (toroidal) velocity dispersion $`\sigma _\phi ^2`$, can be described by the anisotropy distribution

$$a(R,z)1\frac{2\sigma _\phi ^2}{\sigma _\mathrm{m}^2}=\frac{R^2/R_\mathrm{a}^22\alpha }{R^2/R_\mathrm{a}^2+1}.$$
(8)
The velocity dispersion tensor is isotropic when $`a=0`$, tangentially anisotropic when $`a<0`$, and meridionally anisotropic when $`a>0`$. Note that the velocity dispersion anisotropy is constant on cylinders. It can be proved that eq. (8) holds also when the factor $`(QQ_0)^\beta `$ in eq. (1) is replaced by a generic function $`h(QQ_0)`$ (Ciotti et al. 2004). We now impose the system consistency, i.e., we require that

$$\mathrm{}\mathrm{\Psi }=4\pi G\rho ,$$
(9)
with $`\mathrm{\Psi }`$ a nowhere negative function. In order to obtain the numerical solution of eq. (9), the Poisson equation is first recast in dimensionless form, by referring to the physical scales $`R_\mathrm{a}`$ for lengths and $`Q_0`$ for potentials. After rescaling, we are then led to the solution of the following problem:

$$\stackrel{~}{\mathrm{}}\varphi =\lambda \stackrel{~}{\rho },\stackrel{~}{\rho }g(\stackrel{~}{R})(\varphi _\mathrm{T}1)^\gamma \mathrm{\Theta }(\varphi _\mathrm{T}1),$$
(10)
where $`\gamma \alpha +\beta +3/2>0`$, and the dimensionless quantities are defined as $`\varphi _\mathrm{T}\mathrm{\Psi }_\mathrm{T}/Q_0`$, $`\stackrel{~}{R}R/R_\mathrm{a}`$, and $`g(\stackrel{~}{R})\stackrel{~}{R}^{2\alpha }/(1+\stackrel{~}{R}^2)^{\alpha +1/2}`$, with

$$\lambda \frac{16\pi ^2Gf_0R_\mathrm{a}^{2\alpha +2}Q_0^{\alpha +\beta +1/2}B(\alpha +3/2,\beta +1)2^{\alpha +1/2}}{2\alpha +1}.$$
(11)
It follows that

$$\rho =\frac{Q_0}{GR_\mathrm{a}^2}\frac{\lambda \stackrel{~}{\rho }}{4\pi }$$
(12)
and

$$\stackrel{~}{\sigma }_\mathrm{m}^2=\frac{2}{\gamma +1}(\varphi _\mathrm{T}1)\mathrm{\Theta }(\varphi _\mathrm{T}1),\stackrel{~}{\sigma }_\phi ^2=\frac{2\alpha +1}{(\gamma +1)}\frac{(\varphi _\mathrm{T}1)\mathrm{\Theta }(\varphi _\mathrm{T}1)}{1+\stackrel{~}{R}^2}.$$
(13)
Equation (10), a classical non-linear elliptic partial differential equation, is then solved by using a Newton iteration scheme, described in detail by Ciotti et al. (2004). In general, for assigned model parameters we found convergence only for finite intervals of $`\lambda `$. In Figure 1 we show the meridional sections of the density distribution of two representative truncated and fully self-consistent ($`\mathrm{\Psi }_{\mathrm{ext}}=0`$) toroidal models, obtained by fixing $`\alpha >0`$ (see eqs.  and ), for two different values of $`\lambda `$.

The work of G.B. and L.C. has been partially supported by a grant CoFin 2000 of the Italian MIUR.


text:
# Role of structural fluctuations in the insertion into complex host matrices

## I Introduction

Intercalation processes find their application in many technologically important domains, such as the design of hydrogen-storage systems, rechargeable high-energy batteries, electrochromic devices, (see Ref. for a review), electroactive polymers, and superconductors . Microscopically the insertion is a complex process involving many effects (e.g. the charge transfer, elastic response of the matrix, the permselectivity, etc). Even well-characterized materials, like crystalline or layered compounds, exhibit rather complicated elastic properties, involving restructuring, staging and random distortion of the galleries. For amorphous (or porous) matrices  the situation is complicated by the host heterogeneity, that may change in the course of intercalation. For instance, structural heterogeneity of electroactive polymers changes because of the polymer swelling. Complex geometry and poorly characterized energetics make it practically impossible to construct a mechanical model. Therefore, statistical mechanical investigations in this domain face difficulties, induced by the system complexity. In this situation it seems reasonable to develop a new theoretical scheme, that would be capable of making reasonable predictions under restricted information.

For this purpose we combine the standard statistical mechanics and the information theory approach. The latter has been introduced as an alternative (with respect to the standard statistical mechanics) tool for the statistical description of many-body systems. Such an approach has been successfully applied to a description of non-equilibrium steady states and the self-organization in various complex (physical, chemical, biological, etc.) systems. One of the main advantages of the information theory is its capability of describing the equilibrium and nonequilibrium states in the framework of a unified scheme, resembling the usual Gibbs statistics, provided that a suitable information (entropy) measure is chosen. On the other hand, the approach does not rely upon microscopic details. This allows one to study the systems deviating from the standard equilibrium conditions. Quenched systems, like fluids adsorbed in porous media, spin glasses  or liquid glassformers could serve as relevant examples.

In this paper we consider complex systems in which the information at a microscopic level is not sufficient for making a link between the system mechanics (Hamiltonian) and its statistics (probability distribution). Then the we can split the system into subsystems with different levels of description. One is a dynamic subsystem, which evolves according to a Hamiltonian, containing unknown parameters. The latter reflect a coupling to a stochastic subsystem which is governed by a probability distribution. This differs from the usual quenching in, at least, two aspects. The dynamic subsystem can influence the stochastic one and their coupling determines both the thermodynamic behavior and the shape of the probability distribution. For instance, by the analogy with the intercalation system, we can expect a change of the matrix volume upon the fluid absorption. Secondly, one necessarily deals with two entropic impacts: the thermodynamic entropy (due to the dynamic subsystems) and the information entropy (due to the probability distribution). Constructing from these two a suitable entropy measure, we investigate the strain distribution and the guest thermodynamics in homogeneous and heterogeneous matrices.

## II Maximum entropy estimation

Let us assume that the states of a system are labelled by a continuous variable $`x`$, with $`p(x)`$ being the probability of a state $`x`$. Nevertheless, because of the system complexity, the probabilities are not known a priori. The only available information is on a constrained quantity

$$Q=𝑑xp(x)q(x)$$
(1)
which is the expectation value of a physical observable (e.g. energy, heat flux, etc.). In addition the probability distribution is assumed to be normalized

$$𝑑xp(x)=1$$
(2)
Given this information can we estimate another relevant quantity $`G=𝑑xp(x)g(x)`$? For this purpose one has to restore the probability distribution $`p(x)`$.

The basic role in this approach is played by the entropy $`S`$ that gives a measure of missing information concerning the system state. If we have the entropy as a functional of the probability distribution $`S=S[p(x)]`$, then, according to the information theory approach, the probability can be estimated through the entropy maximization under constraints (1),(2). If we take the entropy (or the information) measure in the Shannon form

$$S=𝑑xp(x)\mathrm{ln}p(x)$$
(3)
then such a variation procedure gives

$$p(x)=\frac{e^{\gamma q(x)}}{Z};Z=𝑑xe^{\gamma q(x)}$$
(4)
where $`\gamma `$ is a Lagrange multiplier which has to be determined from the constraint (1). For instance, if $`Q`$ is associated with the equilibrium internal energy, then $`p(x)`$ becomes the conventional Gibbs distribution, with $`\gamma =1/(kT)`$.

## III Complex systems

As is discussed above, by complexity we mean a situation when the available microscopic information on a system is not sufficient for its description in terms of a Hamiltonian. Thus, only one subsystem, say $`\{s\}`$, is governed by a Hamiltonian $`H[\sigma ,s]`$, while the rest (the surrounding) is specified by a set of relevant parameters $`\{\sigma \}`$ which appear with a probability distribution $`P(\sigma )`$. Therefore, we have a coupling of a dynamic system to a stochastic one, which however influence each other. This means that $`P(\sigma )`$ is unknown, the only available information concerns with an observable quantity

$$E=𝑑\sigma P(\sigma )E(\sigma )$$
(5)
that gives, for instance, a characteristic energy scale for the surrounding. In the case of spin glasses the counterparts of $`\{s\}`$ and $`\{\sigma \}`$ systems are the spin variables and the random fields (or exchange constants), respectively. For porous media, $`\{\sigma \}`$ can be identified with the matrix (e.g. pore size) and $`\{s\}`$ \- with the adsorbate degrees of freedom.

For a given configuration $`\{\sigma \}`$ one can calculate the partition function of the dynamic subsystem

$$Z(\sigma )=(ds)e^{\beta H[\sigma ,s]}$$
(6)
where $`\beta =1/(kT)`$ is the inverse temperature. Then all the thermodynamic characteristics are known. For instance, the free energy

$$F(\sigma )=\frac{1}{\beta }\mathrm{ln}Z(\sigma )$$
(7)
and the internal energy

$$U(\sigma )=\frac{d}{d\beta }\mathrm{ln}Z(\sigma )$$
(8)
determine the thermodynamic entropy

$$S_T(\sigma )=\beta [U(\sigma )F(\sigma )]$$
(9)
Then the relevant thermodynamic quantities can be obtained by averaging over all realizations of $`\sigma `$ . For instance, the average free energy is given by

$$F=F(\sigma )=(d\sigma )P(\sigma )F(\sigma )$$
(10)
Note that, in contrast to the conventional equilibrium, now we have two entropies - the thermodynamic entropy

$$S_T=(d\sigma )P(\sigma )S_T(\sigma )$$
(11)
and the one related to the information on the probability distribution

$$S_I=(d\sigma )P(\sigma )\mathrm{ln}P(\sigma )$$
(12)
For a given distribution $`P(\sigma )`$ the value of $`S_T`$ is the average thermodynamic entropy of the coupling $`\{s\}`$-$`\{\sigma \}`$. If, however, the distribution is not known, then $`S_T`$ can be viewed as an additional entropy measure for inferring $`P(\sigma )`$ through the variation procedure.

In the case when the surrounding affects the dynamic subsystem, but not vice versa (pure quenching) $`P(\sigma )`$ can be determined through the $`S_I`$ maximization under the constraint (5) and the normalization condition. This gives

$$P(\sigma )=\frac{e^{\gamma E(\sigma )}}{Z_0}Z_0=𝑑\sigma e^{\gamma E(\sigma )}$$
(13)
where $`E(\sigma )`$ can be estimated through model or even scaling arguments. For instance, $`E(\sigma )=a\sigma ^2`$ is the energy to create a field of the magnitude $`\sigma `$. As is discussed above, the Lagrange multiplier $`\gamma `$ should be determined from the condition (5). Once the probability distribution is determined, the thermodynamic quantities can be calculated straightforwardly.

Such a decoupling of the two subsystems is however an idealization, which could be more or less acceptable depending on a concrete situation. In general, the surrounding should respond to the evolution of the dynamic subsystem. For instance, the porous materials could change the volume upon accommodation of the fluid, like in the case of intercalation systems. Also, it has been demonstrated  that the effect of the quenching on the thermodynamics of $`HCl`$-ice interfaces weakens with increasing ice film thickness.

If the subsystems do influence each other, then the probability distribution must be evaluated starting from the total entropy $`\mathrm{\Sigma }(S_T,S_I)`$, that describes the uncertainty concerning the state of the overall system. Obviously, $`\mathrm{\Sigma }(0,S_I)=S_I`$ and $`\mathrm{\Sigma }(S_T,0)=S_T`$ when there is no loss of information on one of the subsystems. The main problem is to construct $`\mathrm{\Sigma }(S_T,S_I)`$ for nonvanishing $`S_T`$ and $`S_I`$. For instance, we can formally expand around the perfectly ordered state $`S_I=0,S_T=0`$

$$\mathrm{\Sigma }(S_T,S_I)=\mathrm{\Sigma }(0,0)+\kappa S_T+\lambda S_I+\delta S_TS_I+\mathrm{}$$
(14)
As a first approximation the total entropy $`\mathrm{\Sigma }(S_T,S_I)`$ can be estimated as a quasi-additive combination (we drop $`\mathrm{\Sigma }(0,0)`$ as an irrelevant constant)

$$\mathrm{\Sigma }(S_T,S_I)=S_I+\kappa S_T$$
(15)
where $`\kappa `$ is a parameter reflecting the coupling between the dynamic and the stochastic subsystems. Then the distribution function can be inferred by maximizing (15) under the constraint (5). This leads to

$$P(\sigma )=\frac{e^{\gamma ^{}E(\sigma )+\kappa S_T(\sigma )}}{Z}Z=𝑑\sigma e^{\gamma ^{}E(\sigma )+\kappa S_T(\sigma )}$$
(16)
where the Lagrange multiplier $`\gamma ^{}`$ should be determined from the constraint (5). It is known  that the entropy $`S_T(\sigma )`$ determines the probability of thermodynamic fluctuations for a given value of $`\sigma `$ according to

$$P_\phi (\sigma )=\frac{e^{\kappa S_T(\sigma )}}{𝑑\sigma e^{\kappa S_T(\sigma )}}$$
(17)
It should be emphasized that $`P_\phi (\sigma )`$ does not describe the fluctuation of $`\sigma `$, but the fluctuation of thermodynamic variables related to $`\{s\}`$-subsystem (e.g. temperature, density, etc) for a given value of $`\sigma `$. Then the resulting distribution (16) can be represented as

$$P(\sigma )=P_\phi (\sigma )\frac{e^{\gamma ^{}E(\sigma )}}{e^{\gamma ^{}E(\sigma )}_\phi }$$
(18)
where the average $`\mathrm{}_\phi `$ is taken over the thermodynamic fluctuations of the dynamic subsystem. We thus see that the statistics involves two ingredients: the ”energetic”, $`E(\sigma )`$ and the entropic $`S_T(\sigma )`$ or $`P_\phi (\sigma )`$. On the other hand, all thermodynamic quantities (like $`S_T=S_T(\sigma )`$) are also affected by these impacts.

## IV Strain distribution in the course of insertion

In this section we analyze the interplay of these two impacts considering various aspects related to insertion systems. In this case the guest plays the role of $`s`$-subsystem, whose coupling to the host matrix is given by the matrix strain $`\epsilon `$ which is a counterpart of $`\sigma `$-subsystem. We are interested in determining the strain distribution $`P(\epsilon )`$ and the guest thermodynamics $`S_T`$, taking into account the host-guest coupling $`S_T(\epsilon )`$.

### A Homogeneous matrix

As an illustration of our approach we consider the simplest example of guest insertion into a ”homogeneous” host matrix. ”Homogeneous” means that the matrix properties (the structure, the binding site distribution, etc.) are similar for different matrix domains (e.g. graphite galleries). It is well-known that the matrix becomes strained with increasing guest concentration $`x`$. For a given $`x`$ the actual strain is a sum of the domain strains $`_i\epsilon _i`$. For homogeneous matrices we may expect that all the domains are equivalent, such that $`_i\epsilon _i=n\epsilon `$, where $`n`$ is the number of domains.

Therefore, the host-guest coupling leads to a strain distribution $`P(\epsilon )`$ and affects the guest thermodynamics. The strain-dependent guest entropy $`S_T(\epsilon )`$ can be expanded around $`\epsilon =0`$

$$S_T(\epsilon )=S_T(0)+\alpha n\epsilon +\mathrm{}$$
(19)
Here $`S_T(0)`$ is the guest entropy in the case when the matrix is not strained. This corresponds to the lattice gas model describing the guest adsorption on a rigid lattice. Note that all the quantities above ($`S_T(0),\alpha ,\epsilon `$) are concentration dependent. The strain distribution $`P(\epsilon )`$ can be determined by maximizing $`S_I+\kappa S_T`$ under the constraint that the average strain

$$E=n\epsilon =n𝑑\epsilon P(\epsilon )\epsilon $$
(20)
is known (e.g. measured experimentally). It can be estimated as the sample volume dilatation (or $`c`$-axis expansion for layered compounds).

According to eq. (18) we arrive at

$$P(\epsilon )=\frac{e^{(\gamma ^{}\kappa \alpha )\epsilon }}{𝑑\epsilon e^{(\gamma ^{}\kappa \alpha )\epsilon }}$$
(21)
Determining $`\gamma ^{}`$ from the constraint (20) we obtain

$$P(\epsilon )=\frac{e^{\epsilon /E}}{𝑑\epsilon e^{\epsilon /E}}$$
(22)
It should be noted that exactly the same distribution would be recovered if only $`S_I`$ is used as the entropy measure. This is a consequence of what we called homogeneity, that is the average $`E`$ is a good estimation of the actual strain $`\epsilon `$, which almost does not fluctuate. Then the experimentally available information (20) is relevant to the thermodynamic functions of interest. In other words, if the experimental data allow one to recover an ”exhaustive” picture of the system, then any refinement of the entropy measure (like introducing $`S_I+\kappa S_T`$ instead of $`S_I`$) gives no new information.

Based on the probability distribution the guest thermodynamics can be predicted using the experimental data on the concentration dependence of the strain $`E=E(x)`$

$$S_T=S_T(0)+\alpha E$$
(23)
Having determined $`S_T`$, we can recover all the thermodynamic functions in the standard way. The free energy is also linear in the strain (for simplicity we take the same coefficient $`\alpha `$): $`F=F(0)+\alpha E`$. Based on this we obtain the guest chemical potential as

$$\mu =\mu (0)+\alpha \frac{dE}{dx}+E\frac{d\alpha }{dx}$$
(24)
This results coincides with the one obtained in our earlier studies  combining the lattice gas model with the continuum elasticity theory. The latter implies that the stress can be defined as $`\alpha =dF/dE`$. If the stress $`\alpha `$ is composition independent, then from (24) we recover the well-known result of Larche and Cahn .

### B Heterogeneous matrix

In reality the host matrices are not homogeneous. Even for materials which are conventionally regarded as well-characterized, the heterogeneity is practically unavoidable in the experimental conditions (e.g. because of defects, domain structures, etc). On the other hand the amorphous or porous matrices are essentially heterogeneous and their properties are known only statistically. In this case it is difficult to ”invent” a simple experimental test (like (20)), giving an exhaustive description. We assume that the insertion of guest species into a heterogeneous matrix results in a strain distribution for various matrix domains. For simplicity we do not consider spatial strain correlations. This means that the domain strains $`\epsilon _i`$ are randomly distributed according to $`P(\{\epsilon _i\})=_iP(\epsilon _i)`$ independently of the domain position $`i`$.

From the experimental point of view (e.g. measuring the matrix strain $`E`$ or from another suitable test) we know only the range in which the individual $`\epsilon _i`$ can vary. The range of $`\epsilon _i`$ variation should correlate with the average $`E`$. For instance, if $`E=0`$, then $`\epsilon _i`$ are symmetrically distributed around $`0`$. For concreteness we choose: $`0<\epsilon _i<g_i(E)`$. Then, maximizing $`S_I`$ under the normalization condition for $`P(\epsilon _i)`$, we obtain the following step-wise distribution

$$P_0(\epsilon _i)=\frac{H(\epsilon _i)H(g_i\epsilon _i)}{g_i}$$
(25)
where $`H(x)`$ is the Heaviside step function. The distribution $`P_0(\epsilon _i)`$ is the one we can recover using the experimentally available information.

The guest thermodynamics can be predicted as before

$$S_T=S_T(0)+\underset{i}{}\alpha _i\epsilon _i_0$$
(26)
$$\mu =\mu (0)+\underset{i}{}\alpha _i\frac{d}{dx}\epsilon _i_0$$
(27)
where $`\epsilon _i_0=g_i/2`$ and $`\alpha _i`$ is assumed to be composition independent. Note that the results above are based on the estimation of the internal strains $`\epsilon _i_0`$, using the experimentally available information on the ”external” (observable from outside) strain $`E`$. As a simple approximation we may suppose $`g_i(E)=f(S_i)E/n`$, that is, each $`\epsilon _i_0`$ contributes to the observed strain $`E`$ proportionally to the domain size $`S_i`$, with $`_i\epsilon _i_0=E`$.

These results can be refined if we explicitly take into the host-guest coupling and determine $`P(\epsilon _i)`$ using $`S_I+\kappa S_T`$. This leads to the distribution

$$P(\epsilon _i)=\frac{e^{\kappa \alpha _i\epsilon _i}}{\underset{0}{\overset{g_i}{}}𝑑\epsilon _ie^{\kappa \alpha _i\epsilon _i}}$$
(28)
which is more informative in the sense that now not all values of $`\epsilon _i`$ are equally probable. The average

$$\epsilon _i=\frac{1+e^{\kappa \alpha _ig_i}(\kappa \alpha _ig_i1)}{\kappa \alpha _i(e^{\kappa \alpha _ig_i}1)}$$
(29)
is different from $`\epsilon _i_0`$ obtained above. It is instructive to analyze the origin of this difference. For this purpose we expand $`\epsilon _i`$ in terms of the coupling $`\kappa `$

$$\epsilon _i=\epsilon _i_0+\alpha _i\kappa \mathrm{\Delta }_i^0+\mathrm{}$$
(30)
where

$$\mathrm{\Delta }_i^0=\left[\epsilon _i^2_0\epsilon _i_0^2\right]=g_i^2/12$$

describes the strain fluctuations. The latter are related to the domain size fluctuations and thus can be called the structural fluctuations. In the case of homogeneous matrices all the domains are equivalent, then the structural fluctuations vanish $`\mathrm{\Delta }_i^0=0`$ and $`\epsilon _i=\epsilon _i_0`$. Therefore, for heterogeneous matrices our estimation of the strain statistics becomes more sensitive to structural fluctuations if we change the entropy measure. This reflects the fact that for heterogeneous matrices the actual strain $`\overline{E}`$ is different from the one ($`E`$) measured as the sample dilatation. Their difference

$$\overline{E}E=\underset{i}{}[\epsilon _i\epsilon _i_0]=\underset{i}{}\alpha _i\kappa \mathrm{\Delta }_i^0$$
(31)
takes into account the internal deformations which are not detectable experimentally. This result is quite general, it does not rely upon the simple form for $`P_0(\epsilon _i)`$ and $`P(\epsilon _i)`$ discussed here as examples. The only essential fact is that these two are different (compare (13) and (18)). The difference has its physical origin in the mutual coupling of the guest thermodynamics and the host structure. The coupling makes it possible to ”restore” (at least partially) the information inaccessible from experimental tests.

Based on the refined distribution (28) we can calculate the guest thermodynamics.

$$S_T=S_T(0)+\underset{i}{}\alpha _i\epsilon _i_0+\kappa \underset{i}{}\alpha _i^2\mathrm{\Delta }_i^0+\mathrm{}$$
(32)
Considering $`S_T`$ as a measure of the guest thermodynamic fluctuations (see above), we conclude that the latter are coupled to the host structural fluctuations. This reflects the fact that the matrix response to the guest insertion is not only due to the host nature but also involves the host- guest coupling. Comparing the chemical potential

$$\mu =\mu (0)+\underset{i}{}\alpha _i\frac{d}{dx}\epsilon _i_0+\kappa \underset{i}{}\alpha _i^2\frac{d}{dx}\mathrm{\Delta }_i^0$$
(33)
with (27), obtained using the experimental information, we see that (33) gives a more precise estimation, but also contains essential physics. In particular, if $`_i\epsilon _i_00`$, then from (27) one would conclude that the insertion is almost topotactic and that the internal host distortion is irrelevant. In contrast, eq. (33) implies that although the matrix is not strained on average, the expanded and contracted domains do contribute to the guest energetics. This is observed experimentally on various disordered host matrices, which exhibit negligible volume expansion in comparison to well-structured materials. Nevertheless, the electrochemical response of disordered insertion systems is remarkably different. This can be understood in terms of two effects. First is the configurational transitions described by $`\mu (0)`$. Structural disorder induces a broad distribution of binding energies on different matrix sites. This leads to quite steep insertion isotherms with multiple short plateaus. A simple model capable of recovering this effect can be found in ref. .

The second effect is due to the structural transitions, which are supplemented by structural fluctuations (see eq. (33)). In order to analyze this contribution within our simple model, let us recall that $`\epsilon _i_0=g_i(E(x))/2`$ and $`\mathrm{\Delta }_i^0=g_i^2(E(x))/12`$. Assuming the simplest non-trivial concentration dependence $`g_i(E(x))=b_ix`$, we obtain from (33):

$$\mu =\mu (0)+\delta \mu +Kx$$
(34)
where $`\delta \mu =_ib_i\alpha _i/2`$ induces a shift of the isotherm in comparison to the strain free state, and $`K=\kappa _i\alpha _i^2b_i^2/6`$ plays the role of additional repulsive ($`K>0`$) interaction. Since the chemical diffusion coefficient is proportional to $`d\mu /dx`$, it is clear that the structural fluctuations tend to increase diffusivity in comparison to the strain free state. This seems to be a general tendency for disordered matrices, as it follows from experimental studies.

## V Conclusion

The coupling of structural and thermodynamic fluctuations in the course of various-type insertion processes is investigated within a combination of Gibbsian statistics and the information theory approach. The main focus is on complex (heterogeneous) insertion systems, where the microscopic information is not sufficient for making a link between the system mechanics (Hamiltonian) and its statistics (probability distribution). Then the system can be viewed as an entropic coupling of a dynamic and a stochastic subsystem. It is shown that the coupling makes it possible to restore (at least partially) the information, inaccessible from experimental tests. This enables one to make physically reasonable predictions under limited information on the system.

It should be emphasized that our results are also applicable to matrices which are initially homogeneous but become heterogeneous in the course of insertion. For instance, this takes place in $`Li`$-graphite intercalation compounds, which exhibit staging phenomena. In this case the occupied and empty graphite galleries are not equivalent. Also, if a gallery is partially occupied, then one deals with distorted domains. In this case the strain fluctuations as well as their spatial correlation should be taken into account. In a simple approximation this physics is reflected by our result (see eq. (33)). Our approach is also capable of describing the swelling of electroactive polymers in the course of intercalation. We plan to address these issues in a future study.