---
abstract: 'A simple model for the angular dependent interaction between $C_{60}$ molecules in face centered cubic lattice is proposed and analyzed by use of the rigorous bifurcation approach.  The quantitative results for the orientational phase transition and the characteristics of the ordered phase are in good agreement with the experimental data.' 
address: ' Institute for High Pressure Physics, Russian Academy of Sciences, Troitsk 142092, Moscow region, Russia'
author:
- 'T.I.Schelkacheva and E.E.Tareyeva'
title: ' Orientational Phase Transition in Solid $C_{60}$: A Model Approach'
---

The orientational phase transition in solid $C_{60}$ is of much current interest.  The orientational ordering has been a subject of extensive experimental investigations (see, e.g.,  [@dav; @hein; @dav1; @blinc; @blas; @braj] ); some theoretical researches were performed, too  [@mich; @yha; @lu; @heid; @lap; @laun; @gam].  However, [*ab initio*]{} calculations fail to reproduce the experimental results. 

In this Letter we develop a simple model for the angular dependence of the intermolecular potential in solid $C_{60}$. The model is based on the ideas of preferred orientations due to David et al.   [@dav] and to Lapinskas et al.   [@lap], and on the maximal exploit of symmetry considerations.  We apply to this model interaction the rigorous approach based on the Lyapunov–Schmidt theory of bifurcation of solutions of nonlinear integral equations and obtain quantitative results for the transition temperature and the distribution of molecular orientations in the ordered phase.  These results occur to reproduce the experimental data. 

As is established in a number of experiments $C_{60}$ crystallizes in a face centered cubic (fcc) structure.  At ambient temperature the molecules rotate almost freely with centers on the fcc lattice sites, so that the space group is $Fm{\bar 3}m$ (see, e.g.,  [@hein]).  When the temperature decreases to $T_S \approx 260 K$ the first order orientational phase transition takes place: the sites of the initial fcc latice become divided between four simple cubic sublattices (see fig.1)

(180,140)(0,20) (90,95)[(0,2)[30]{}]{} (90,95)[(2,0)[30]{}]{} (90,95)[(-1,-1)[22]{}]{} (90,95) (90,135)[$z$]{} (125,95)[$y$]{} (65,65)[$x$]{} (100,100)[0]{} (60,60)[(0,2)[37]{}]{} (60,113)[(0,2)[37]{}]{} (60,105) (57,103)[A]{} (30,30)[(0,2)[37]{}]{} (30,83)[(0,2)[37]{}]{} (27,73)[A]{} (30,75) (150,60)[(0,2)[37]{}]{} (150,113)[(0,2)[37]{}]{} (150,105) (147,103)[A]{} (120,30)[(0,2)[37]{}]{} (120,83)[(0,2)[37]{}]{} (120,75) (117,73)[A]{}

(60,60)[(2,0)[37]{}]{} (113,60)[(2,0)[37]{}]{} (105,60) (102,58)[B]{} (30,30)[(2,0)[37]{}]{} (83,30)[(2,0)[37]{}]{} (75,30) (72,28)[B]{} (30,120)[(2,0)[37]{}]{} (83,120)[(2,0)[37]{}]{} (75,120) (72,118)[B]{} (60,150)[(2,0)[37]{}]{} (113,150)[(2,0)[37]{}]{} (105,150) (102,148)[B]{}

(150,60)[(-1,-1)[10]{}]{} (130,40)[(-1,-1)[10]{}]{} (135,45) (132,43)[D]{} (60,60)[(-1,-1)[10]{}]{} (40,40)[(-1,-1)[10]{}]{} (45,45) (42,43)[D]{}

(60,150)[(-1,-1)[10]{}]{} (40,130)[(-1,-1)[10]{}]{} (45,135) (42,133)[D]{}

(150,150)[(-1,-1)[10]{}]{} (130,130)[(-1,-1)[10]{}]{} (135,135) (132,133)[D]{} (70,5)[Fig.1.]{} 

with own preferable molecular orientation in each sublattice.  The broken symmetry space group is $Pa{\bar 3}$.

Moreover, the neutron–diffraction experiments  [@dav] have shown that the orientations in the ordered state are so that the electron–rich regions (the interpentagon double bonds) face the electron–deficient regions of the neighboring $C_{60}$ molecule: the centers of pentagons or the centers of hexagons.  It was shown  [@dav; @dav1; @blas] that the ratio of the number of molecules in those two states is about $60:40$ at the phase transition temperature and increases when the temperature decreases.  This remaining orientational disorder is usually believed to cause the orientational glass transition at $T_G \approx 90K$ now confirmed by various experimental technics (see, e.g.,  [@glass]).  These two minima of the intermolecular angle dependent energy were obtained by numerical calculations and were shown to be much lower than the energies of other mutual orientations of the pair of molecules (see, e.g.,  [@yha; @lu; @gam; @yha1]).  In those calculations the previously obtained charge distribution for the isolated $C_{60}$ molecule  [@mol] was taken into account.  Usually recent calculations use the intermolecular potential of Sprik et al.   [@sprik]: a sum of 6-12 and Coulomb interactions between 60 atoms $C$ and 30 double–bond centers $D$ and between each other: $$\begin{aligned}
\Phi(1,2) & = &
\sum_{k\in C(1)} \sum_{k'\in C(2)}
4 \epsilon \left \{ \left( \frac {\sigma_{CC}} {R_{kk'}}
\right)^{12} - \left( \frac {\sigma_{CC}}
{R_{kk'}}\right)^6\right\} \nonumber\\
& + &\sum_{k \ne k', k, k' \in C,D}
4 \epsilon \left \{ \left( \frac {\sigma_{CD}} {R_{kk'}}
\right)^{12} - \left( \frac {\sigma_{CD}}
{R_{kk'}}\right)^6\right\} \nonumber\\
& + &\sum_{k\in D(1)} \sum_{k'\in D(2)}
4 \epsilon \left \{ \left( \frac {\sigma_{DD}} {R_{kk'}}
\right)^{12} - \left( \frac {\sigma_{DD}}
{R_{kk'}}\right)^6\right\} \nonumber\\
& + & \sum_{k,k'\in C,D} \frac {q_{k} q_{k'}} {R_{kk'}}
\label{full}\end{aligned}$$ Here $\epsilon = 1.293 meV, \sigma_{CC} = 3.4 \AA, \sigma_{CD} = 3.5
\AA, \sigma_{DD} = 3.6 \AA, q_D = - 0.35e, q_C = - q_D/2.$

Rigorously speaking we are interested in the angular part of this complicated interaction represented in terms of multipole–multipole interaction of point–like multipoles on the sites of rigid fcc lattice with coefficients to be calculated from  (\[full\]).  The general form of this angular part is $$\Phi_{ij} (\omega_i, \omega_j) =
\sum_{l;\nu, \tau}
C_{\nu, \tau}^l (\omega_{ij})  u _{l \nu} (\omega_i) u_{l \tau}
(\omega_j),
\label{mult}$$ with $l = 6, 10, 12, 16, 18, ...$ due to the icosahedral molecular symmetry $I_h$. In the Eq.  (\[mult\]) $\omega_i$ are the angles describing the orientation of the molecule on site $i$, for example, Euler angles and $u_{l \nu }$ – some kind of harmonics.  However, we simplify the problem and develop a model orientational interaction.  As the angular dependent interaction is rather short–ranged we can restrict ourselves by the nearest–neighbor interactions. 

We follow the main ideas of the papers  [@dav; @lap] and use the restricted number of allowed orientations instead of free continuous rotations.  Let us take into account in the energy  (\[mult\]) only the orientations with pentagons, hexagons or double bonds directed towards 12 nearest neighbors in fcc lattice.  The $C_{60}$ molecule is constructed in such a way that if 6 of its 12 pentagons (or 6 of its 20 hexagons) face 6 nearest neighbors double bonds ($P$ and $H$ states of Lapinskas et al.  [@lap]), then 6 of its 30 interpentagon double bonds face the remaining 6 nearest neighbors.  Now the energy matrix elements can take only three values: $J_0$ – the energy of the general mutual position, $J_P$ – pentagon [*versus*]{} double bond and $J_H$ – hexagon [*versus*]{} double bond.  These energies in our model can be compared with those calculated in  [@yha; @lu; @gam; @yha1] as functions of the angular displacements of the molecule at (0,0,0).  Following  [@lap], and putting $J_0 = 0$ we obtain from the fig.2(b) of the paper  [@yha] $J_P = - 300K$ and $J_H = - 110K$. Now we leave the paper  [@lap] and follow our own way. 

The energy matrix elements $J_P$ and $J_H$ connect the states of molecules only in the allowed orientations.  So, only allowed linear combinations of $u_{l \nu}$ enter the Eq.   (\[mult\]).  The theoretical curve in  [@yha] makes no difference between the number $l$ of harmonics and describes the effect of all of them.  So, in the framework of our model calculation it is possible to build up the allowed functions using only the harmonics with $l=6$: we need only their transformation properties.  We restrict ourselves to $l=6$, however the coefficients $J_P$ and $J_H$ are not some of $C_{\nu ,\tau }^6$ given in  (\[mult\]) but effectively take into account higher order terms.  Let us construct the functions $P_i(\omega )$ and $H_i(\omega )$ explicitly in terms of cubic harmonics $K_{m} \equiv K_{6,m}, m = 1, 2, ..., 13$ (see, e.g.,  [@heid]).  All functions $P_i$ and $H_i$ are the sums of $K_m$, invariant under the icosahedral symmetry of the molecule (i.e. belonging to the $A_{1g}$ representation of the icosahedral group $I_h$) if icosahedrons are naturally oriented in one of 8 properly chosen coordinate systems.  The states $P_i (H_i)$ have 6 pentagons (hexagons) and 6 double bonds directed towards 12 nearest neighbors along different $[100]$ axes.  $P_1(\omega )$ describes the molecule rotated from the standard orientation $B$ (following  [@harr]) about $[111]$ axis through the angle $97.76125^o$. The angle for $H_1(\omega )$ is $37.76125^o$. The functions $P_2 (\omega
)$,$P_3 (\omega )$ and $P_4 (\omega )$ (or $H_2,
H_3, H_4$) are obtained from $P_1(\omega )$ $(H_1)$ by subsequent counter–clockwise rotations of the molecule by $90^o$ around $z$ axes. 

If written in the standard coordinate frame with Cartesian axes along the cube sides these functions have the following explicit form: $$\begin{aligned}
P_{1}(\omega ) & = &
\alpha_{P} K_{1} (\omega ) + \beta_{P} [K_{8}(\omega ) + K_{9}
(\omega ) + K_{10} (\omega )] \nonumber \\ & + & \gamma_{P}
[K_{11}(\omega ) + K_{12} (\omega ) + K_{13} (\omega )],\nonumber\\
P_{2}(\omega ) & = & \alpha_{P} K_{1} (\omega ) +
\beta_{P} [ - K_{8}(\omega ) + K_{9} (\omega ) - K_{10} (\omega
)] \nonumber\\& + & \gamma_{P} [ - K_{11}(\omega ) + K_{12}
(\omega ) - K_{13} (\omega )], \nonumber\\ P_{3}(\omega ) & = &
\alpha_{P} K_{1} (\omega ) + \beta _{P} [K_{8}(\omega ) - K_{9}
(\omega ) - K_{10} (\omega )] \nonumber\\& + & \gamma _{P}
[K_{11}(\omega ) - K_{12} (\omega ) - K_{13} (\omega )],
\nonumber\\
P_{4}(\omega ) & = & \alpha_{P} K_{1} (\omega ) + \beta _{P}
[ - K_{8}(\omega ) - K_{9} (\omega ) + K_{10} (\omega )]
\nonumber\\& + &
\gamma _{P}
[ - K_{11}(\omega ) - K_{12} (\omega ) + K_{13} (\omega )],
\label{ph}\end{aligned}$$ with $\alpha_{P}=-0.38866;$ $\beta _{P}=0.31486;$ $\gamma _{P}=-0.42877.$ The functions $H_i(\omega )$ have the same form as $P_i(\omega )$ but with the coefficients $\alpha _{H}= 0.46588;$ $\beta _{H}=0.37740;$ $\gamma _{H}=0.34432$. The functions are normalized to unity. 

Let us now treat our model by use of bifurcation approach in the mean–field approximation.  As is well known the mean–field approach often brings one to the formulation of the broken space symmetry problem in terms of the bifurcation of solutions of nonlinear integral equations for distribution functions (see, e.g., the review  [@koz]).  In particular, the bifurcation approach was used in the case of orientational phase transitions in molecular crystals in the Refs.  etc.  The simplified version was originally developed by James and Keenan for solid methane  [@jkee] and by Michel, Copley and Neumann  [@mich] and by Heid  [@heid] for solid $C_{60}$. We shall follow our papers on hydrogen  [@tt1; @ttdan; @tt2; @sch]. 

In the mean–field approximation from the first equation of BBGKY hierarchy for the orientational distribution functions or by minimizing the orientational free energy one can obtain the following nonlinear integral equation [@ttdan]: $$g_i ( \omega_i) + \frac{1} {\Theta} \sum_{i \ne j} G_j
\int d \omega_j \Phi_{ij} ( \omega_i, \omega_j)
 e^{g_i (\omega_j)} = 0;
\label{main}$$ $g_i (\omega_i) = \ln[\frac {f_i (\omega_i)} {G_i}]$, $f_i ( \omega_i)$ – one–particle orientational distribution function for a molecule on $i$-th lattice site, the constants $
G_i$ are the normalization constants. 

In our case of solid $C_{60}$ where there are four sublattices (see fig.1) and four kinds of unknown distribution functions we obtain from  (\[main\]) the following system of four nonlinear integral equations:

$$\begin{aligned}
g_1(\omega )& + &\lambda \int d \omega ' [
B(\omega ,\omega ') G_2 e^{g_2(\omega ')} + A(\omega ,\omega ')
G_3 e^{g_3 (\omega ')}
\nonumber\\& + &
D(\omega ,\omega ') G_4 e^{g_4 (\omega
')}] = 0,
\nonumber\\
g_2(\omega )& + &\lambda
\int d \omega ' [ B(\omega ,\omega ') G_1 e^{ g_1(\omega ')} +
A(\omega ,\omega ') G_4 e^{g_4 (\omega ')}
\nonumber\\& + &
D(\omega ,\omega ')
G_3 e^{g_3 (\omega ')}] = 0,\nonumber\\
g_3(\omega )& + &\lambda \int d \omega ' [ B(\omega ,\omega ')
G_4 e^{g_4(\omega ')} + A(\omega ,\omega ') G_1 e^{g_1 (\omega
')}
\nonumber\\& + &
D(\omega ,\omega ') G_2 e^{g_2 (\omega ')}] = 0,
\nonumber\\
g_4(\omega )& + &\lambda \int d \omega'
[ B(\omega ,\omega ') G_3 e^{g_3(\omega ')} + A(\omega ,\omega ')
G_2 e^{g_2 (\omega ')}
\nonumber\\& + &
D(\omega ,\omega ') G_1 e^{g_1 (\omega
')}] = 0. 
 \label{int}\end{aligned}$$

Here $\lambda=1/T$, $A(\omega,\omega')$, $B(\omega ,\omega ')$, $D(\omega ,\omega ')$ are the sums of interactions over nearest neighbors in the sublattices $A$, $B$ and $D$ (see fig.1), respectively.  For example, the sum in the plain perpendicular to the $x$ axis can be written explicitly in the form: $$\begin{aligned}
D(\omega ,\omega ') & = & \nonumber\\
=2 \{[(P_1 (\omega ) &+& P_4 (\omega )]
J_{P} + (H_1 (\omega ) + H_4 (\omega _)] J_H]\nonumber\\
 \times  [
P_2 (\omega ')& +& P_3 (\omega ') + H_2 (\omega ') + H_3 (\omega
')]\nonumber\\
+  [P_2 (\omega ) &+& P_3 (\omega ) + H_2
(\omega ) + H_3(\omega )] \nonumber\\
\times  [(P_1 (\omega ') &+&
P_4 (\omega ')) J_P + (H_1 (\omega ') + H_4 (\omega ')) J_H]
\nonumber\\
 +  [(P_2 (\omega ) &+& P_3 (\omega )) J_P + (H_2
(\omega ) + H_3 (\omega )) J_H]\nonumber\\
\times  [ P_1 (\omega ') &+& P_4 (\omega ') + H_1 (\omega ') +
H_4 (\omega ')] \nonumber\\
 +  [P_1(\omega ) &+& P_4 (\omega )
+ H_1 (\omega ) + H_4 (\omega )]\nonumber\\
 \times
[(P_2(\omega ') &+& P_3(\omega ')) J_P + (H_2 (\omega ') + H_3
(\omega ')) J_H]\}
\label{pph}\end{aligned}$$ and analogously for two other sublattices. 

The equations  (\[int\]) are well known Hammerstein equations  [@ham].  In the case of finite domain of integration when the fixed point principle is valid there exists detailed theory for such equations (see  [@kra]).  We use the standard methods (see, e.g.  [@tren]).  At high temperature the system  (\[int\]) has only trivial solution $g_i( \omega_i)=0$, corresponding to the orientationally disordered phase.  At the bifurcation points $\lambda _{\alpha}$ new solutions with broken symmetry appear ( $\lambda _{\alpha} > 0$).  For $\lambda = \lambda _ {\alpha} (1+ \mu)$ the functions $g_i^{\alpha }( \omega_i)$ can be written as series in integer or fractional powers of $\mu$. These powers are defined by the bifurcation equation (see  [@tren]) corresponding to the system  (\[int\]).  In our case we have $$g_i (\omega ) = \mu h_i (\omega ) + \mu^2 x_i (\omega ) +
...$$ because among the integrals $\int d \omega K_{m_1} K_{m_2}
K_{m_3}$ with $m=8,...,13$ there are some which are not equal to zero.  This means the first order phase transition  [@tt1; @ttdan; @tt2]. 

The bifurcation points are the eigenvalues $\lambda _ \alpha$ of the linearized system corresponding to  (\[int\]): $$\begin{aligned}
h_1(\omega ) &+& \frac {\lambda} {4 \pi} \int d \omega ' [
B(\omega ,\omega ') h_2(\omega ') + A(\omega ,\omega ') h_3
(\omega ') \nonumber\\
& + &D(\omega ,\omega ') h_4 (\omega ')] = 0.\nonumber\\
\dots &\dots & \hbox to 6cm{\dotfill}
\label{sys}\end{aligned}$$ The functions $h_i$ can be written in the form $$h_i(\omega) = \sum_{\nu} h_{i}^{\nu} K_{\nu}
(\omega), \label{hk}$$ so that the eigenvalues $\lambda _ \alpha $ define not only the bifurcation temperatures but the relations between nonzero coefficients $h_{i}^{\nu}$ (that is the symmetry of the new phase), too. 

In the case of the full interaction one can obtain all possible broken symmetry phases compatible with the initial symmetry and the condition of positive temperature value (see, e.g., the case of hydrogen  [@tt1; @tt2; @sch]).  Now we have truncated the interaction and reduced the problem.  Nevertheless, there remain still two quantitative characteristics we aim to obtain: the bifurcation temperature and the relation between the weights of $P$ and $H$ functions in the solution. 

Using  (\[hk\]) it is easy to rewrite the system  (\[sys\]) as the system of linear algebraic equations for the coefficients $h_i^{\nu}$ with $i=1,2,3,4$ and $\nu = 1,8,9,...,13$.

Using the explicit form of the matrices $A, B, D$ it is easy to obtain the only nonzero elements: $$A_{1,1}  =  B_{1,1} = D_{1,1} \equiv u,$$ $$A_{8,8}  =  B_{9,9} = D_{10,10} \equiv v,$$ $$A_{11,11}  =  B_{12,12} = D_{13,13} \equiv z,$$ $$A_{8,11} =  A_{11,8} = B_{9,12} = B_{12,9} = D_{10,13} =
D_{13,10} \equiv w.$$ One can write the elements $u, v, z$ and $w$ in terms of the coefficients $\alpha _P,  \beta
_P, \gamma _P, \alpha _H, \beta _H, \gamma _H$ and energies $J_P, J_H$ and obtain the following values: $u = 32 \cdot
5.046, v = 32 \cdot 94.127, z = 32 \cdot 7.665, w = - 32
\cdot 37.155$.

The determinant of the algebraic system is factorized in $2 \times 2$ determinants, so that the eigenvalues $\lambda _ \alpha $ can be easily obtained.  Among the values $\lambda _ \alpha $ there are two positive values.  The first one $\lambda_1  = 4 \pi /u$ corresponds to the solution proportional to $K_1$ and is of no interest now.  The second is the positive solution of the equation $$1 - \frac {\lambda} {4 \pi} (v + z) + \frac
{{\lambda}^2} {(4 \pi)^2} (vz - w^2) = 0,
\label{eql}$$ namely $ \lambda_b = 0.00364 K^{-1}$ or $ T_b = 275 K$. The corresponding nontrivial eigenfunctions have $h_i^1 = 0,$ and the other coefficients $h_i^{\mu }$ are subject to some constraints.  If we add the condition for the functions $h_i(\omega )$ to transform one into another under the action of the cubic group rotation elements which leave the fcc lattice invariant, then only three of the coefficients remain to be independent and the functions $h_i$ can be written in the following form: $$\begin{aligned}
h_1(\omega ) & = & a P_1 (\omega ) + b H_1 (\omega ) + c K_1
(\omega),\nonumber\\
h_2(\omega ) & = & a P_3 (\omega ) + b H_3 (\omega ) + c K_1
(\omega),\nonumber\\
h_3(\omega ) & = & a P_4 (\omega ) + b H_4 (\omega ) + c K_1
(\omega),\nonumber\\
h_4(\omega ) & = & a P_2 (\omega ) + b H_2 (\omega ) + c K_1
(\omega),
\label{hpk}\end{aligned}$$ $$a \alpha _P + b \alpha _H + c = 0,
\label{fin1}$$ $$a \beta _P + b \beta _H = Q (a \gamma _P + b \gamma _H),
\label{fin2}$$ $$Q = \frac {1 - v \frac
{\lambda} {4 \pi}} { \frac {\lambda} {4 \pi}w}.$$

Using the numerical value for $Q$ we obtain immediately: $$\rho _P =\frac {a} {a + b} = 0.608; \qquad \rho _H=\frac {b} {a
+b} = 0.392,
\label{ro}$$ so that not only the transition temperature, but also the ratio of the number of molecules in $P$ and $H$ states occur to coincide with the experimental data  [@dav; @dav1; @blas] $\rho _P:\rho _H=60:40.$

To obtain the remaining unknown coefficient we use the equations of the second order in $\mu$ ($ g_i (\omega ) = \mu h_i (\omega ) + \mu ^2 x_i (\omega );
\lambda  = \lambda _c (1 + \mu )$).  The system has the form: $ \hat L \times {\bf x} = {\bf R} $, where $\hat L$ is the linear $4 \times 4$ operator $\hat L \times {\bf h} = 0.$ The system has nontrivial solutions for $x_i$ if the right hand side is orthogonal to the solutions for $h_i$ obtained before.  All 16 equations $$\int d\omega R_i (\omega ) h_j (\omega )
= 0
\label{fin3}$$ are identical due to the symmetry of coefficients.  Solving equations  (\[fin1\]),  (\[fin2\]) and  (\[fin3\]), we obtain finally $ a = - 25.7;\quad b =  - 16.6;\quad c = -2.27.$ The minus sign means that the solution goes in the direction of higher temperatures ($ \mu =
- \tau, \tau = (T - T_b)/T_b$).  The solution has the turning point $T_t$ which is some Kelvins higher than $T_b$. The actual first order phase transition obtained from the free energy behavior takes place between these two points.  The details of this calculation are to be published elsewhere. 

To conclude, we developed a simple model for angle dependent interaction for $C_{60}$ molecules in the fcc cubic lattice.  We used rigorous analytic approach based on the Lyapunov–Schmidt theory of bifurcation of solutions of nonlinear integral equations to treat this model.  As result we obtained the first order phase transition, the bifurcation temperature $ T_b = 275 K$, the $Pa{\bar 3}$ symmetry of the ordered phase and the ratio $\rho $ of the number of molecules with pentagon facing neighbor double bond to the number of molecules with hexagon facing neighbor double bond  (\[ro\]) near the phase transition in good agreement with the experimental data. 

This work was partially supported by Russian Foundation for Basic Researches (Grant No.  98-02-16805). 

The authors would like to thank V.N.Ryzhov and V.A.Davydov for useful discussions. 

W.I.F.David, R.M.Ibberson, T.J.S.Dennis, J.P.Hare, and K.Prassides, Europhys.  Lett., [** 18**]{}, 219 (1992).  P.A.Heiney,G.B.M.Vaughan,J.E.Fischer,N.Coustel,\
D.E.Cox,J.R.D.Copley,D.A.Neumann,\
W.A.Kamitakahara,K.M.Creegan,D.M.Cox,\
J.P.McCauley, Jr.,A.B.Smith III, Phys. Rev., [**B 45**]{}, 4544 (1992). 

W.I.F.David, Europhys. News, [**24**]{}, 71 (1993).  R.Blinc, J.Selinger, J.Dolinsek, and D.Arcon, Phys. Rev., [**B 49**]{}, 4993 (1994). 

O.Blaschk, G.Krexner, Ch. Maier,R.Karawatzky,\
Phys. Rev. [**B 56**]{}, 2288 (1997). 

V.V.Brazh kin, and A.G.Liapin, UFN (Uspechi, Moscow) [**166**]{}, 893 (1996).  K.H.Michel, J.R.D.Copley, and D.A.Neumann,\
Phys. Rev.Lett. [** 68**]{}, 2929 (1992).  T.Yildirim and A.B.Harris, Phys. Rev. [ **B 46**]{}, 7878 (1992).  J.P.Lu, X.-P.Li, and R.M.Martin, Phys. Rev.Lett., [** 68**]{}, 1551 (1992).  R.Heid, Phys. Rev. [**B 47**]{}, 15912 (1993).  S.Lapinskas, E.E.Tornau, and A.Rosengren, Phys. Rev.[**B 49**]{}, 9372 (1994).  P.Launois, S.Ravy, and R.Moret, Phys. Rev.[**B 55**]{}, 2651 (1997).  Z.Gamba, Phys. Rev. [**B 57**]{}, 1402 (1998).  F.Gugenberger,R.Heid,C.Meingast,P.Adelman,M.Braun,\
H.Wühl,M.Haluska,and H.Kuzmany, Phys. Rev.Lett. [** 69**]{}, 3774 (1992).  T.Yildrim, A.B.Harris, S.C.Erwin, and M.R.Pederson,\
Phys. Rev. [**B 48**]{}, 1888 (1993).  B.P.Feuston, W.Andreoni, M.Parinello,\
and E.Clementi, Phys. Rev. [**B 44**]{}, 4056 (1991);\
N.Troullier and J.L.Martins, Phys. Rev. [**B 46**]{}, 1754 (1992), and references therein.  M.Sprik, A.Cheng, and M.L.Klein, J.Phys. Chem. [** 96**]{}, 2027 (1992).  A.B.Harris, Physica [**A 205**]{}, 154 (1994).  J.J.Kozak, Adv. Chem. Phys., [** 40**]{}, 229 (1979).  E.E.Tareyeva, and T.I.Trapezina, Physics Letters, [**A51**]{}, 114 (1975); [*ibid*]{} , [**A60**]{}, 217 (1977).  E.E.Tareyeva, and T.I.Trapezina, DAN USSR (Soviet Phys.  Doklady), [**223**]{}, 823 (1975); Theor. Math. Phys. ( Moscow), [**26**]{}, 269 (1976).  E.E.Tareyeva and T.I.Schelkacheva, Theor. Math. Phys. ( Moscow), [**31**]{}, 359 (1977); V.N.Ryzhov, E.E.Tareyeva and T.I.Schelkacheva, [*ibid* ]{} [**40**]{}, 269 (1979).  H.L.Lemberg, and S.A.Rice, Physica, [**63**]{}, 48 (1973).  A.N.Shramkov, and P.P.Shtifanuk, Theor. Math. Phys. ( Moscow), [**90**]{}, 246 (1992).  T.I.Shchelkacheva, Physics Letters, [**A239**]{}, 397 (1998).  H.M.James,and T.A.Keenan,J.Chem. Phys.[**31**]{},12(1959).  A.Hammerstein, Acta Math., [** 54**]{}, 117 (1930).  M.A.Krasnosel’skii, Topological methods in the theory of Nonlinear Integral Equations, Pergammon press, London, 1964.  M.Vainberg and V.A.Trenogin, Russian Math.  Surveys, [**17**]{}, 1 (1962); M.Vainberg, [*Variational Methods for the Study of Nonlinear Operators*]{}, Holden Day, San Francisco, 1964.