format: text

text:
# Entropy and Temperature of a Quantum Carnot Engine

## 1 Introduction

A classical thermodynamic heat engine converts heat energy into mechanical work by using a classical-mechanical system in which a gas expands and pushes a piston in a cylinder. Such a heat engine obtains its energy from a high-temperature heat reservoir. Some of the energy taken from this reservoir is converted to mechanical work. A heat engine is not perfectly efficient, so some of the energy taken from the heat reservoir is not converted to mechanical energy, but rather is transferred to a low-temperature reservoir (Planck 1927).

A classical heat engine running between a high-temperature reservoir and a low-temperature reservoir achieves maximum efficiency if it is reversible. While it is impossible to construct a working heat engine that is perfectly reversible, in the early 19th century Carnot proposed a mathematical model of an ideal heat engine that is not only reversible but also cyclic (Carnot 1824). The Carnot engine consists of a cylinder of ideal gas that is alternately placed in thermal contact with high-temperature and low-temperature heat reservoirs whose temperatures are $`T_H`$ and $`T_C`$, respectively.

Instead of a classical system of gas, we consider here a quantum-mechanical system of consisting a single particle in contact with a reservoir. The system is described by a time-dependent cyclic Hamiltonian. The statistical ensemble of many such identically-prepared systems that we consider here is characterised by a density matrix. We assume that the system interacts weakly with its environment. It is known that it is possible to extract work from such a system (e.g., Geva $`\&`$ Kosloff 1994; Kosloff, et al. 2000).

In particular, if the evolution of the density matrix is reversible, then we can construct a quantum-mechanical Carnot engine (Bender, et al. 2000). The purpose of the present paper is to investigate the properties of quantum Carnot engine. This is of importance because it improves our understanding of the relationship between thermodynamics and quantum mechanics, an area of considerable interest in quantum theory (e.g., Leff $`\&`$ Rex 1990).

This paper is organised as follows: First, we explain the concept of a quantum engine by generalising the two-state model considered by Bender, et al. (2000) to an infinite-state square-well model and derive equations of states for isoenergetic and adiabatic processes. These results lead naturally to the quantum analogue of the Clausius equality for a reversible cycle. We show that, unlike the result in classical thermodynamics, the Clausius relation obtained here is not based on the change of entropy.

We then introduce the von Neumann entropy and obtain the maximum entropy state subject to isoenergetic requirements. We demonstrate that the maximum von Neumann entropy is consistent with the thermodynamic definition of entropy, as well as the equations of states for quantum Carnot cycle. As a consequence, we are able to determine the temperature of the energy bath directly from two of the diagonal components of the density matrix. We conclude by discussing several open problems.

## 2 Quantum Carnot cycle

We can construct a simple quantum engine using a single particle of mass $`m`$ confined to an infinite one-dimensional square-well potential whose volume (width) is $`V`$. For any fixed $`V`$ it is easy to solve the time independent Schrödinger equation to determine the energy spectrum of the system:

$`E_n(V)={\displaystyle \frac{\pi ^2\mathrm{}^2n^2}{2mV^2}}.`$ (1)
We assume that the width of the well initially is $`V=V_1`$ and that the initial energy of the system is a fixed constant $`E_H`$. The initial state $`\psi (x)`$ of the system is a linear combination $`\psi (x)=a_n\varphi _n(x)`$ of the energy eigenstates $`\varphi _n(x)`$. Thus,

$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}p_nE_n(V_1)=E_H,`$ (2)
where $`p_n=|a_n|^2`$, and $`E_H`$ is bounded below by $`E_HE_1(V_1)`$. Note that the pure state $`\psi (x)`$ characterises a typical element of the ensemble, whose statistical property is thus determined by the density matrix, given, in the energy eigenstates, by $`\rho _{mn}=p_n\delta _{mn}`$.

Starting with the above initial configuration, we define a quantum cycle as follows: First, the well expands isoenergetically; that is, the width of the well increases infinitely slowly while the system is kept in contact with an energy bath. Note that the quantum adiabatic theorem (Born $`\&`$ Fock 1928) states that if the system were isolated during such an expansion, the system would remain in its initial state. That is, the absolute values of the expansion coefficients $`|a_n|`$ would remain constant. Thus, if the system were isolated, the energy of the system (the expectation value of the Hamiltonian) would decrease like $`V^2`$.

However, during this expansion, we simultaneously pump energy into the system in order to compensate this decrease of the energy. Thus, the well, which can be viewed as a one-dimensional piston, expands in such a way that the expectation value of the Hamiltonian is held constant by exciting higher energy states. This implies that, unlike the isothermal process in classical thermodynamics, the temperature of the ‘bath’ in an isoenergetic process is in fact changing. The explicit volume dependence of temperature is derived in Section 5.

During such an isoenergetic expansion, mechanical work is done by the force (one-dimensional pressure) $`P`$ on the walls of the well. The contribution to this force from the $`n`$th energy eigenstate is

$`f_n={\displaystyle \frac{\pi ^2\mathrm{}^2n^2}{mV^3}}.`$ (3)
Hence the force $`P`$ is given by the expectation value $`P=p_nf_n`$. Using the relation $`f_n=2E_n/V`$ and Eq. (2), we obtain the equation of state during an isoenergetic process:

$`PV=2E_H.`$ (4)
This result is identical to the corresponding equation of state for an isothermal process of a classical ideal gas, if we make the identification $`2E_HkT_H`$. Note that the expansion coefficients $`a_n`$ of the wave function change as a function of the width $`V`$ while the well expands isoenergetically from $`V_1`$ to $`V_2`$.

Second, the system expands adiabatically. During an adiabatic process the eigenstates $`\varphi _n(x)`$ change as a function of $`V`$, but the values $`|a_n|`$ remain constant. Therefore, the expectation value of the Hamiltonian $`E=p_nE_n(V)`$ decreases during the process because each $`E_n`$ decreases with increasing $`V`$ while all $`p_n`$ are kept fixed. The force $`P`$ in this case is determined by differentiating the energy $`E`$ with respect to $`V`$. Thus, the equation of state during an adiabatic process for the square-well potential is

$`PV^3=2V_2^2E_H,`$ (5)
which is a quantum analogue of the corresponding equation for a classical ideal gas. The system expands adiabatically until its volume reaches $`V=V_3`$. At this point the expectation of the Hamiltonian decreases to $`E_C`$. Because the squared coefficients $`p_n`$ of the wave function remain constant during an adiabatic process, the value of $`E_C`$ is given by $`E_C=(V_2^2/V_3^2)E_H`$.

Following the adiabatic expansion, the system is compressed isoenergetically until $`V=V_4`$, with the expectation value of the Hamiltonian fixed at $`E_C`$, and finally it is compressed adiabatically until the width of the system returns to its initial value $`V=V_1`$. This cycle is reversible, and we find that the efficiency of the quantum engine is given by

$`\eta =1{\displaystyle \frac{E_C}{E_H}},`$ (6)
a formula analogous to the classical thermodynamic result $`\eta =1T_C/T_H`$ of Carnot (1824).

## 3 Quantum Clausius relations

We have demonstrated above an example of a quantum system with a slowly-changing time-dependent Hamiltonian from which work can be extracted. The key concept introduced here is an energy bath that maintains the expectation value of the Hamiltonian. This idea can be developed further to establish the quantum counterparts of the classical thermodynamic relations. To begin with, let us first consider the theorem of Clausius.

During an isoenergetic expansion, the amount of energy transferred to the system to maintain the expectation of the Hamiltonian is determined by the integral

$`Q_H={\displaystyle _{V_1}^{V_2}}𝑑VP(V)=2E_H\mathrm{ln}(V_2/V_1),`$ (7)
where $`P(V)=2E_HV^1`$. The amount of energy absorbed during an isoenergetic compression can be determined similarly with the result $`Q_C=2E_C\mathrm{ln}(V_3/V_4)`$. Thus, for a closed, reversible cycle we obtain the quantum Clausius equality

$`{\displaystyle \frac{Q_H}{E_H}}+{\displaystyle \frac{Q_C}{E_C}}=0,`$ (8)
because $`V_2/V_1=V_3/V_4`$ for a closed cycle.

In classical thermodynamics the Clausius equality states that the total change $`𝑑S`$ of entropy in a reversible cycle is zero, where the differential of entropy $`dS=dQ/T`$ is the ratio of the absorbed heat and the bath temperature. Therefore, the relation (8) suggests that in quantum mechanics the entropy change in an isoenergetic process might be given by $`dQ/E`$, the ratio of the absorbed energy to the bath energy. However, as we show below, this quantity does not determine the amount of entropy change, and hence the quantum Clausius relation is not a condition for entropy. Instead, it merely implies that the total energy absorbed, for given bath energies, must vanish in a reversible cycle. If the cycle is irreversible, then we have $`𝑑Q/E<0`$.

## 4 Entropy maximisation

To understand the change of entropy in a quantum Carnot cycle we consider the von Neumann entropy

$`S={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}p_n\mathrm{ln}p_n`$ (9)
associated with the density matrix $`\rho _{mn}`$, expressed in terms of the energy eigenstates. Let us first consider the isoenergetic process and introduce a scaling parameter $`\lambda `$ such that $`V=\lambda V_1`$. Then, during an isoenergetic process the probability $`p_n`$ that the system is in the $`n`$th eigenstate changes in $`\lambda `$, subject to the constraints

$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}p_n=1\mathrm{and}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}n^2p_n=\lambda ^2,`$ (10)
where we have chosen units in which $`\pi ^2\mathrm{}^2/(2m)=1`$ and taken the initial condition $`E_H=E_1(V_1)`$. There are infinitely many combinations of $`p_n`$ satisfying these constraints, each of which determines the entropy. Hence, quantum mechanically it appears that any one of such states is allowed in an isoenergetic process. However, there is a unique density matrix that satisfies the thermodynamic requirement of the change of entropy, and this is the state that maximises the von Neumann entropy.

To show this, we must determine the density matrix that maximises the entropy subject to the constraints (10) and then obtain the associated entropy $`S`$. We can perform this maximisation directly without using Lagrange multipliers. The two constraints in (10) imply that two of the diagonal components, say $`p_k`$ and $`p_l`$, of the density matrix are determined. Therefore, if we differentiate the constraints (10) with respect to $`p_n`$, we find that

$`{\displaystyle \frac{p_k}{p_n}}+{\displaystyle \frac{p_l}{p_n}}=1\mathrm{and}k^2{\displaystyle \frac{p_k}{p_n}}+l^2{\displaystyle \frac{p_l}{p_n}}=n^2.`$ (11)
Solving these linear equations, we obtain

$`{\displaystyle \frac{p_k}{p_n}}={\displaystyle \frac{n^2l^2}{k^2l^2}}\mathrm{and}{\displaystyle \frac{p_l}{p_n}}={\displaystyle \frac{n^2k^2}{l^2k^2}}`$ (12)
for all $`k,ln`$. Substituting these expressions into $`S/p_n=0`$ and solving for $`p_n`$, we get

$`p_n=p_l\left({\displaystyle \frac{p_k}{p_l}}\right)^{(n^2l^2)/(k^2l^2)},`$ (13)
which is also valid for any choice of $`kl`$.

Given the maximum entropy state (13), we can now reexpress the constraints in the form

$`\lambda ^2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\alpha ^{n^2}={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}n^2\alpha ^{n^2},`$ (14)
where we have defined $`\alpha =(p_k/p_l)^{1/(k^2l^2)}`$. This is obtained by substituting (13) in (10), solving for $`p_l`$, and equating the resulting expressions. Similarly, using the parameter $`\alpha `$, the state represented in (13) can be expressed as

$`p_n={\displaystyle \frac{\alpha ^{n^2}}{_{m=1}^{\mathrm{}}\alpha ^{m^2}}}.`$ (15)
Therefore, the corresponding von Neumann entropy is

$`S=(l^2\lambda ^2)\mathrm{ln}\alpha \mathrm{ln}p_l.`$ (16)
This is independent of $`l`$ because of the identity

$`p_k\alpha ^{k^2}=p_l\alpha ^{l^2}.`$ (17)
Thus, given the width scale $`\lambda =V/V_1`$ of the potential well, the relation (14) determines the value of $`\alpha `$, from which we can determine both $`p_n`$ and $`S`$ as functions of the single length scale parameter $`\lambda `$.

We can now analyse the entropy associated with the isoenergetic process. We note first that the change of entropy associated with (16) does not agree with the quantity $`dQ/E`$. More specifically, we verify that only asymptotically for large values of $`\lambda `$, that is, large energies, the two quantities are proportional to each other. In particular, because the value of $`S`$ in (16) is independent of the choices of $`k`$ and $`l`$, we can set $`k=2`$ and $`l=1`$ without loss of generality. Then, asymptotically we have $`p_1p_2`$ for large $`\lambda `$. Furthermore, the maximum entropy state $`\{p_n\}`$ satisfies $`p_k>p_l`$ for all $`k<l`$. Therefore, $`\alpha 1`$ and for large values of $`\lambda `$ we use the asymptotic relation

$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1ϵ)^{n^2}{\displaystyle \frac{\sqrt{\pi }}{2\sqrt{ϵ}}}{\displaystyle \frac{1}{2}}(ϵ0^+)`$ (18)
to determine the entropy change. The result gives $`\mathrm{\Delta }S\mathrm{ln}\lambda `$ when the width changes from $`V_1`$ to $`V=\lambda V_1`$. On the other hand, from (7) we deduce that $`dQ/E=2\mathrm{ln}\lambda `$ for any value of $`\lambda 1`$. This establishes the asymptotic proportionality. However, $`dQ/E`$ does not in general determine the entropy change.

## 5 Equilibrium-state temperature

Next, we show that the maximum von Neumann entropy (16) is consistent with thermodynamic considerations. Recall that an isoenergetic process must proceed infinitely slowly in order that the energy bath at any point during the process be in thermal equilibrium. Thus, this process requires an infinite amount of time to be realised. Then from the principles of statistical mechanics (Schrödinger 1952) we deduce that, if we set $`V_1=1`$, the temperature of the bath is given by $`T=1/(\lambda ^2\mathrm{ln}\alpha )`$. This is because we can make the identification $`\alpha ^{n^2}=\mathrm{exp}(E_n/T)`$ in Eq. (15). On the other hand, temperature is given in thermodynamics by the formula

$`{\displaystyle \frac{1}{T}}={\displaystyle \frac{dS}{dQ}},`$ (19)
where $`dQ=PdV`$. In the present consideration we have $`dQ=2\lambda ^1d\lambda `$ from the equation of state $`PV=2E_H`$ and the relation $`dQ=PdV`$, while from (16) we find that $`dS/d\lambda =2\lambda \mathrm{ln}\alpha `$. Therefore, the thermodynamic definition of temperature in (19) gives

$`T={\displaystyle \frac{1}{\lambda ^2\mathrm{ln}\alpha }},`$ (20)
which agrees with the statistical mechanical consideration above. We conclude therefore that the maximum von Neumann entropy determines the density matrix of the system in an isoenergetic process. Furthermore, the temperature of the energy bath can be obtained from (20), where $`\alpha `$ is determined by any two diagonal components $`p_k`$, $`p_l`$ of the density matrix. In particular, the temperature of the bath changes continuously in $`\lambda `$ during an isoenergetic process, even though the energy of the system is held fixed.

We now consider the adiabatic process. In this case, there is no net energy absorbed by the system, so the thermodynamic entropy remains constant. On the other hand, because the diagonal components $`p_n`$ of the density matrix are constant during an adiabatic process, the von Neumann entropy (9) also remains constant, in agreement with the thermodynamic consideration.

These results can be expressed by the phase diagram in the entropy-volume plane, as shown in figure 1. The diagram shows that the thermodynamic constraint $`V_2/V_1=V_3/V_4`$ is satisfied by, and only by, the maximum von Neumann entropy. Thus, we emphasise that the quantum state corresponding to the maximum von Neumann entropy, as obtained by microscopic analysis, is entirely consistent with the thermodynamic equations of states.

## 6 Discussion

In the analysis above we have considered only the square-well potential model, which can be interpreted as the quantum analogue of the classical ideal gas. A natural extension of this work would be to consider other Hamiltonians. For example, for a harmonic potential whose energy eigenvalues are $`E_n=(n+\frac{1}{2})\mathrm{}\omega `$, a dimensional argument shows that the characteristic length scale is given by $`1/\sqrt{\omega }`$. Therefore, the equations of states for the harmonic potential become identical to those for the square-well potential, although the volume dependence of the entropy is different from (16). For other Hamiltonians, however, the equations of states are not in general identical to those obtained here, and it would be interesting to find explicit examples of other models exhibiting nonideal behaviours.

Another important issue is to understand whether there is any role played by the geometric phases, a concept that does not have an analogue in the classical Carnot cycle. Although the quantum Carnot engine is indeed cyclic in terms of $`p_n`$, the coefficients $`a_n`$ of the wave function of any one of the representatives in the ensemble inevitably pick up geometric phases as the system goes through a cycle. Therefore, strictly speaking, the wave function is not cyclic in the quantum Carnot cycle. The question then is whether there is any physically observable evidence of the geometric phase in the present context, and if so, how do we represent that in the mixed state context (cf. Uhlmann 1996).

Finally, another interesting idea that should be studied is the possibility of constructing a cyclic quantum engine (or a quantum refrigerator by reversing the cycle) that requires only a finite amount of time to complete. Note that each cycle of a classical Carnot engine is infinitely long. However, in quantum mechanics, it is known, for example, that by choosing a specific form of time dependent Hamiltonians, it is possible to construct a quantum-mechanical cycle in an essentially arbitrary short time scale (Mielnik 1986). It would be of interest to determine if such an idea can be applied to create a finite-time Carnot engine, or whether such dynamics would be incompatible with the requirement of maximising the von Neumann entropy. An analysis on the performance of finite-time cycle in the presence of friction has been studied in Feldmann $`\&`$ Kosloff (2000).

We wish to express our gratitude to A. Howie and M. F. Parry for stimulating discussions. DCB gratefully acknowledges financial support from The Royal Society. This work was supported in part by the U.S. Department of Energy.


text:
# Two-Stage Rotational Disordering of a Molecular Crystal Surface: C60

\[

## Abstract

We propose a two-stage mechanism for the rotational surface disordering phase transition of a molecular crystal, as realized in C60 fullerite. Our study, based on Monte Carlo simulations, uncovers the existence of a new intermediate regime, between a low temperature ordered $`(2\times 2)`$ state, and a high temperature $`(1\times 1)`$ disordered phase. In the intermediate regime there is partial disorder, strongest for a subset of particularly frustrated surface molecules. These concepts and calculations provide a coherent understanding of experimental observations, with possible extension to other molecular crystal surfaces.

\]

Surface-specific thermal disordering phenomena such as for instance surface melting , are of considerable interest. Conceptually, they demonstrate how a chunk of the metastable phase (e.g, the liquid phase within the solid domain of stability) can be imprisoned in full equilibrium by restricting it to an interface region of space. Practically, they provide at a temperature below the triple point a germ of the disordered phase that makes a major difference to all surface-related properties, as well as to the bulk transformation kinetics.

In molecular crystals, thermal surface disordering is richer than in atomic solids. Because of the presence of rotational degrees of freedom, the bulk develops rotational and positional melting transitions that are separate, and both of which might possess a surface counterpart. One is thus led to seek the existence of a well defined surface rotational melting for a molecular crystal. However until recently, very little phenomenology has gathered addressing that issue.

That situation recently changed after the discovery of the C60 molecule and its subsequent availability in macroscopic quantities in its crystalline form, the fullerite. Here the near sphericity of the molecules favors rotational disordering at a relatively low temperature, when positional disorder and evaporation are still irrelevant owing to a sufficiently large cohesion . On the theoretical side, the remarkable symmetry of the molecule leads to considerable simplifications with the help of symmetry adapted functions , such as those originally used for the CD4 molecules .

Bulk measurements showed that C60-fullerite undergoes a first-order rotational melting transition at 260K, characterized by fixed orientations of the molecules in the space group $`Pa\overline{3}`$ at low temperature transforming to a plastic phase with rotational isotropy, space group $`Fm\overline{3}m`$, at high temperature . In the $`Pa\overline{3}`$ phase, the fullerite structure is composed of four sublattices. Neutron powder diffraction studies  of this phase have revealed the orientational configuration of the C60 molecules. The optimal orientation is such that for all adjacent couples of C60, an electron-rich double bond of one molecule faces an electron-poor pentagon of the other. In the $`Fm\overline{3}m`$ phase, the molecules are still not completely free to rotate as they are still correlated, and moreover they feel a crystal field of cubic symmetry. However, in this phase angular correlations are short-ranged, and molecules are all equivalent, as opposed to the $`Pa\overline{3}`$ phase, with four inequivalent molecules and long-range angular correlations.

The structure of a free (111) surface of fullerite can differ from that of a bulk (111) plane due to the lowering of the number of neighbors. The first experiments showing a specific surface disordering were high-resolution reflection electron-energy-loss spectra (EELS) obtained by Goldoni et al . This technique is sensitive to the surface layer only due to the short escape depth of the scattered electrons. The authors reported an abrupt jump around 230K in the energy of the C60 triplet exciton. This discontinuity is an indication that the fullerite (111) surface undergoes the rotational disordering phase transition at a temperature some 30K lower than the bulk. A similar transition temperature ($`230\pm 20`$K) has been observed in LEED studies of C60 films  but the incertitude in the temperature was too large to interpret it as a specific surface effect.

At first sight, such a lowered surface transition temperature is reminiscent of ordinary, positional surface melting of solids, where surface premelting precedes and initiates the bulk transition . However, it soon appeared that in fullerite the surface disordering took place in the topmost layer only , in contrast with standard surface melting, where the liquid film grows progressively invading the crystal. At the same time, a conceptual difference was also noted. Positional order in a crystal is forced exclusively by interparticle interaction; rotational ordering instead involves both intermolecular angular interactions and a cubic crystal field that would be present and effective on each molecule even if all the surrounding molecules were perfect spheres. A surface molecule thus feels not only less rotational interactions, but also a crystal field which may generally frustrate bulk order. On the other hand, due to the weak short range nature of intermolecular potentials, the molecules in the second layer, immediately below the surface are already very close to a bulk situation. The system, between 230K and 260K, can thus be regarded as an ordered solid with a disordered surface.

We now underline the fact that experimental data suggest in fact not one, but two surface transitions, the second one taking place at even lower temperature. There is a less-marked discontinuity around 160K in the triplet exciton energy revealed by EELS . It was explained as a consequence of the molecular rotation relaxation time which at this temperature becomes similar to the exciton lifetime. However an independent evidence for a transition or crossover at that lower temperature is the presence around 160K of a shoulder in the LEED intensity of the $`Pa\overline{3}`$ diffraction spots versus temperature, absent in the bulk X-ray intensity . Helium atom scattering (HAS) diffraction studies of the (111) surface of fullerite by Glebov et al clearly showed a change from a (2X2) to a (1X1) structure taking place at 235K. But, in addition, the behavior of the Debye-Waller Factor (DWF) of the (0,0) Bragg spot between 140K and 235K was completely unexpected. Instead of decreasing monotonically with increasing temperature, the DWF started increasing around 140K and reached a maximum at 235K. This was explained as a consequence of preempted disordering of C60 molecules at low coordination defect sites. This behavior, found to be reversible, cannot be explained if the surface layer were totally ordered below 235K.

In the present work, we studied the surface orientational ordering of fullerite by classical Monte-Carlo (MC) simulation. We considered a C60 molecule as a rigid body, and used the potential developed by Lamoen  to model the intermolecular interactions. That model considers 210 interacting centers on each fullerene: 60 are localized on the carbon atoms, 60 are at the center of the single bonds, and three are distributed along each of the 30 double bonds to take their spatial extension into account. The repulsive part of the intermolecular potential is a sum of Born-Mayer terms between all the interacting centers (c) and the attractive part is made of $`R^6`$ van der Waals interactions between the atoms (a)

$$V=\underset{a_1a_2}{}\frac{A}{r_{a_1a_2}^6}+\underset{c_1c_2}{}B_{c_1c_2}\mathrm{exp}\left[C_{c_1c_2}r_{c_1c_2}\right]$$
(1)
Our MC simulations were based on a standard Metropolis, acceptance-rejection algorithm. The molecules were placed in a semi-infinite fcc lattice with a (111) surface. We considered twelve non-equivalent molecules, four in the first layer (the surface layer), four in the second layer, and four in the third layer. All other molecules in the system were assumed to be copies of this twelve-molecules basic cell, with periodic boundary conditions applied in two dimensions. In the direction of the bulk, the second layer was reproduced at the fourth layer, shifted horizontally in order to simulate the fcc stacking. The center of mass of the molecules were held fixed to lattice sites, while the temperature dependent lattice parameter of the structure was set to the experimental value. The twelve molecular orientations were allowed to change randomly. We performed a complete calculation of the energy at each MC step.

Earlier theoretical studies of the rotational bulk phase transition have shown the usefulness of a set of order parameters based on the cubic harmonics . We used an efficient and original algorithm to generate these cubic harmonics. Starting from generating functions of the irreducible representations $`k`$ ($`A_1`$, $`A_2`$, $`E`$, $`T_1`$, $`T_2`$) of the cubic group $`O_h`$, it is possible to expand them into Taylor series and to extract from that expansion the number of irreducible representations for each manifold $`l`$ of the intermolecular potential. Using the generators $`\mathrm{\Gamma }\left(k\right)`$ of each irreducible representation, one can construct automatically the cubic harmonics

$$\left[\mathrm{\Gamma }^{\left(k\right)}\left(g\right)\right]E^{(l,k)}\left(\stackrel{}{r}\right)=E^{(l,k)}\left(\mathrm{\Gamma }^{\left(T_1\right)}\left(g\right)\stackrel{}{r}\right)$$
(2)
where $`g`$ belong to $`O_h`$ and $`E^{(k,l)}`$ is the integrity basis of order $`l`$ for the irreducible representation $`k`$. Since the cubic harmonics must be homogeneous polynomials of the coordinates, this system can be solved easily, and one obtains analytical expressions for all the cubic harmonics, to be used with the atomic coordinate of the 60 atoms of the molecule. In particular, the following expression was calculated

$$h_n^{(l,k)}=\underset{a=1}{\overset{60}{}}E^{(l,k)}\left(\stackrel{}{r_a^n}\right)$$
(3)
where $`\stackrel{}{r_a^n}`$ is the position vector of atom $`a`$ in the $`n^{th}`$ molecule in the unit cell. Using the C60 atomic coordinates automatically introduces the icosahedral symmetry of the molecules. Another advantage of this method is to naturally allow for a future inclusion of internal molecular vibrations, not treated at this stage. Due to the icosahedral symmetry, $`h_n^{(l,k)}`$ differ from zero for $`l`$ = 6, 10, 12, 16, 18 $`\mathrm{}`$ . We used $`l=10`$ throughout. For each molecule (index $`n`$), we calculated a scalar order parameter

$$\eta _n^{(l,k)}=|h_n^{(l,k)}|^2$$
(4)
where $`h_n^{(l,k)}`$ is an average over the configurations generated by the MC algorithm. For all the irreducible representations of the $`O_h`$ group, this parameter was traced as a function of the temperature to observe the phase transition as it has the largest amplitude in the ordered phase.

As mentioned in  $`T_{1g}`$ and $`T_{2g}`$ are the two primary order parameters. We adopted the irreducible representation $`T_{2g}`$ as the best indicator of the rotational phase transition, as the corresponding order parameter vanishes in the disordered phase . Figure 1 shows the temperature dependence of that order parameter for the twelve molecules of our system. The eight bulk-like molecules of layers 2 and 3 (represented by triangles), plus all their infinite copies, undergo a clear phase transition at 250K. This value is affected by size, in particular by artificial correlations built in a small system like the one we considered. By comparison, the mean-field theory transition temperature was however even smaller, namely 200K . However, only the dominant part of the multiple order parameters was originally used in the mean-field molecular equations of ref. . An extension of the number of order parameters did raise the mean field transition temperature to 230K .

At low temperature, the molecules in the (111) surface (represented by circles in figure 1) keep the same orientation as in the bulk, corresponding to the $`(2\times 2)`$ superstructure observed in HAS diffraction. In this configuration, shown in figure 2, one can immediately distinguish two kinds of molecular orientations. Three molecules in the surface $`(2\times 2)`$ unit cell have only a two-fold axis normal to the surface and the $`Pa\overline{3}`$ threefold axis 30 degrees off. The fourth molecule, which we denote as BS, has instead the three-fold axis normal to the surface. Calculations performed by rotating one molecule while keeping a fixed orientation of all the other molecules show that the “crystal field” angle-dependent potential felt by the rotating molecule is different, and shallower for the BS molecule than for the three others. This suggests that the rotational thermal behavior of the BS molecule must differ from that of the other surface molecules .

In fact with increasing temperature, the four surface molecules do behave differently. The BS practically becomes a quasi-free rotator above around 150K whereas the other three molecules do so only above a much higher phase transition temperature, which is around 230K.

In the low temperature phase, the surface symmetry is threefold due to the local threefold axis of the BS molecule. At 150K, we observe a dramatic reduction of threefold order with strong rotational disordering of the BS molecule. Although strictly speaking the threefold order of that molecule cannot vanish, it does however become so small to assimilate it to a free rotator, in which case the surface symmetry would become sixfold because the BS molecule environment is characterized by a sixfold axis. Thus 150K represents a crossover temperature, where there is no strict change of surface symmetry and probably no proper surface phase transition, but there is nonetheless a change from essential order to essential disorder for BS molecules only.

Hence our system exhibits altogether four phases or regimes (see table I):

(a) a fully ordered $`Pa\overline{3}`$ state (effective $`p3`$ surface symmetry) below 150K. (b) an intermediate regime between 150K and 230K with one essentially disordered and three ordered surface molecules per unit cell ($`p6`$ effective symmetry of the surface layer) over a $`Pa\overline{3}`$ ordered bulk. (c) a disordered surface layer on an ordered $`Pa\overline{3}`$ bulk between 230K and 250K. (d) a fully rotationally disordered $`Fm\overline{3}m`$ ($`p6m`$ symmetry of the surface layer) above 250K.

In the simulations, full thermal disordering of the molecules takes place at the topmost layer only in agreement with evidence from EELS . In addition, the main surface transition temperature (230K) agrees very well with the observations by LEED, EELS and HAS. The intermediate BS-disordered regime between 150K and 230K can explain the weaker experimental anomalies mentioned above, in particular the feature at 160K in the triplet exciton energy, signaling in our interpretation an approximate change of surface symmetry from $`p3`$ to $`p6`$.

We note that in LEED and HAS diffraction, this change is not picked up since the unit cell remains $`(2\times 2)`$ up to 230K. We suggest instead that STM experiments be carried out in the appropriate temperature interval to test in more detail our predicted two-stage surface disordering. Our model provides a plausible explanation for the behavior of the (0,0) Debye-Waller factor in HAS diffraction , its increase between 140 and 235K due to the BS molecule rotating freely in that temperature interval. Due to the resulting partial disordering of the $`(2\times 2)`$ surface cell, the intensity of some diffracted beams is reduced and transferred to the specular spot , which compensates the thermal attenuation of the elastic reflection . Our proposed selective thermal disordering of BS molecules may also provide an explanation of the anomalous features observed near 160K in X-ray powder diffraction intensities. Such anomalies were also observed by David  in high-resolution neutron scattering. If these anomalies were not reported later on, it is probably due to the increasing quality (and then the grain size) of the C60 powder, decreasing the surface/bulk ratio.

In summary, the study of fullerite may prove of paradigmatic importance to a much broader class of molecular crystals. Our study indicates that rotational thermal surface melting of a molecular crystal differs from ordinary, positional surface melting, in one important aspect, with several consequences. The rotational ordering forces comprise both crystal fields and intermolecular interactions, while the former is missing in the positional case. At the surface, both interactions and crystal fields are different from the bulk. When, as in fullerite, the crystal field is quantitatively very important in the bulk, its heavy change at the surface may lead to an early thermal disordering taking place strictly in the first layer, as in all other layers the crystal field is bulk-like. Moreover, crystal fields may differ so much between one surface molecule and another, to cause surface thermal disordering to take place in separate stages, the weaker molecules disordering much earlier than the stronger ones.

The authors acknowledge support through the EU project FULPROP, contract number ERBFMRXCT 970155. This research was also funded by NSF and the PAI P4/10 on reduced dimensionality systems. Work at SISSA was sponsored by MURST COFIN99, and by INFM/F . The authors thanks also P. Senet, V. Meunier, L Henrard, P. Rudolf, S. Modesti, and A. Goldoni for many interesting discussions.