invited papers
Pathintegral approach to anharmonic vibration of solids and solid interfaces
^{a}Department of Chemistry, Graduate School of Science, The University of Tokyo, 731 Hongo, Bunkyoku, Tokyo 1130033, Japan
^{*}Correspondence email: toshi@chem.s.utokyo.ac.jp
The temperature dependence of
(extended Xray absorption fine structure) cumulants was investigated for bulk Cu and a thin film of Cu by means of the pathintegral effective classical potential method. By using the semiempirical embeddedatom method as a potential, agreement between the experiments and calculations is found to be excellent for bulk Cu. In the thin Cu(111) film, anisotropic anharmonic vibration was clearly observed; the outofplane vibration is much more enhanced and more anharmonic than the lateral vibration. The results are semiquantitatively consistent with the previous experimental data on Cu(111)/graphite. Such a vibrational enhancement should be the driving force for the roughening transition and/or the surface premelting at higher temperature. The practical usefulness of the pathintegral effective classical potential method combined with the embeddedatom method is demonstrated.Keywords: EXAFS cumulants; path integral; effective classical potential; lattice dynamics; anharmonicity; thin films; anisotropic vibration.
1. Introduction
) and the perturbation formulae are useful to describe molecules and clusters (Yokoyama et al., 1996, 1997; Yokoyama, 1999). It is, however, rather complicated to apply the perturbation theory to largedimensional systems, such as solids and solid surfaces, especially for higherorder cumulants that describe anharmonicity.
spectroscopy has been applied widely to the determination of local structures of crystalline and noncrystalline solids, amorphous materials and liquids, solutions, catalysts, biological and environmental materials, surface and interfaces, and so forth. is commonly recognized to provide information on coordination numbers and interatomic distances. However, since contains information on radial distribution functions around Xrayabsorbing atoms, one can obtain higherorder information on thermal fluctuation, including anharmonicity. Stimulated by interesting experimental observations by means of temperaturedependent existing theories have been exploited in recent years. The quantummechanical statistical perturbation theory gives analytical formulae of the cumulants that are experimentally obtainable (Fujikawa & Miyanaga, 1993Another theoretical approach is the Feynman pathintegral theory, which treats the motion of atoms in real space like a classical image for particles. The realspace representation should be more suitable to describe the radial distribution function, especially of solids observed by et al., 1992, 1995; Kleinert, 1995). The first application of the PIECP theory to cumulants was conducted for a twobody potential (Fujikawa et al., 1997) and comparisons with experimental data of a diatomic molecule of Br_{2} and of solids such as facecentred cubic (f.c.c.) Kr and Ni, were subsequently performed (Yokoyama, 1998).
when compared with the perturbation theory, which is based on the eigenstate (or reciprocal) representation. Moreover, the pathintegral theory can handle greater anharmonicity than the perturbation theory. Although the pathintegral Monte Carlo method is established, the calculations usually require heavy computational loads. A more practical method is to employ the pathintegral effective classical potential (PIECP) theory (CuccoliIn the former part of this article, a brief review of the PIECP theory is at first presented, including a practical usage of the potential function. The PIECP theory is based on the twobody potential, which is sufficiently accurate to describe raregas solids but is insufficient for metals or semiconductors. It is, however, found that one can apply the embeddedatom method (EAM) (Foiles, 1985) to the PIECP theory without modification of the formulae. The results for bulk Cu calculated using the PIECP theory combined with EAM are subsequently given.
In the latter part of this paper, a surface anisotropic and anharmonic vibration is discussed. An enhancement of vibrational amplitudes and anharmonicity of surface atoms is believed to be the driving force of surface roughening transitions, subsequent surface premelting and consequent bulk melting, the transition mechanisms of which have been a longstanding problem from a microscopic point of view. et al., 2000). In this work, the PIECP calculations of thin Cu(111) films are compared with the experimental observations to demonstrate the practical usefulness of the PIECP theory for temperaturedependent studies of real systems.
provides valuable information on local disorder which cannot be obtained by diffraction techniques. In order to observe a surface anisotropic and anharmonic vibration, we have investigated thin Ni(111) and Cu(111) films grown epitaxically on graphite and have found that the outofplane surface vibration is significantly enhanced and becomes more anharmonic compared with the inplane one (Kiguchi2. Theory
2.1. Path integral effective classical potential method
In this section, the PIECP theory is reviewed briefly. Details of the theory can be found in the literature (Cuccoli et al., 1992, 1995; Kleinert, 1995). According to the Feynman theory, the density function ρ(X) (X is the 3Ndimensional Cartesian coordinate) is expressed as a functional form:
where β = (k_{B}T)^{−1}, is the Planck constant divided by 2π, is the Hamiltonian of the system, Z is the and [X(u)] is the Euclidean action,
In the PIECP approximation, the Euclidean action is solved variationally. For a trial function of [X(u)], a solution of a harmonic oscillator should be the best candidate since thermal properties of solids are being investigated:
where is the average path,
and the force constant matrix F and the scalar potential w are variational parameters.
The Cartesian coordinate X is transformed into a normal coordinate Q through Q = ^{t}UM^{1/2}(X − ), where U is the eigenvector matrix of M^{−1/2}FM^{1/2}. One can thus rewrite
where is the diagonal eigenvalue matrix. The density matrix ρ_{0} of the harmonic oscillator is obtained as
where
α_{k} is the difference between the quantummechanical and classical fluctuations of the phonon mode k. The expected value of a certain physical quantity is given as
where is the effective classical potential,
In order to optimize the variational parameters of w and F, the Jensen–Feynman inequality can be used:
where F and F_{0} are respectively the true and trial free energies. The optimized results are as follows:
and
where denotes the dyadic, the ij component of which is defined as
2.2. Lowcoupling approximation in twobody potential solids
The 3Ndimensional integrals in equations (6) and (8) cannot, however, be numerically calculated in the case of a largedimensional system. We will here employ a further approximation: the socalled lowcoupling approximation. This assumes that w and F are independent of . Moreover, we will for simplicity confine ourselves in a monatomic with an m in which only twobody forces are active. The dynamical matrix D is defined as
where F_{oj} is the 3 × 3 matrix of F for atoms o and j, and R_{oj} is the location vector of atom j with respect to atom o. When the eigenvalues and eigenvectors of the matrix D are written as and (μ = 1, 2 and 3 are the phonon branches), the effective classical potential is obtained as
where is the bare classical potential, u(R_{ij}) is the twobody potential between atoms i and j with the distance R_{ij}, and R_{ij}^{0} is the and are respectively expressed as
and
where is the unit vector of . is the difference of the
Debye–Waller factors between the quantummechanical and the classical schemes. The formula can be adopted for a monatomic with a twobody potential, such as raregas solids.2.3. Application of the embeddedatom method to PIECP
Metallic bonds are regarded as a manybody force between valence electron densities and ion cores, and twobody interaction potentials are inadequate to describe thermal properties of metals. The EAM has extensively been employed in classical simulations of metals, which assumes semiempirical parameters based on the densityfunctional theory. The EAM potential V is written as
where ρ_{h,i} is the host electron density,
and is the electron density of atom j at the site of atom i. is the embedding energy of atom i into the host and corresponds to an attractive manybody force between valence electrons and ion cores. is in principle a functional but is replaced by a normal function F_{i} under the localdensity approximation. is the repulsive potential between ion cores i and j. can be calculated by using atomic wavefunctions, and the function forms of F_{i} and φ_{ij} are empirically parameterized.
The EAM potential seems not to be applicable to the PIECP formula (15) since it contains the manybody force of F_{i}. When the EAM potential is Taylor expanded up to the second order, however, we obtain
where is given in the equilibrium geometry and
The expansion includes a very important consequence: the EAM potential is two body within the harmonic approximation. We can thus employ equation (15) of a twobody potential system within the EAM. This is ascribed to the absence of angulardependent terms in the EAM potential, although it is inherently a manybody potential. We can replace the bare classical potential in equation (15) with the EAM potential V in equation (18) for PIECP calculations of metals.
3. Computational details
For the PIECP calculations of bulk Cu, the above formalism is available without further approximation. The employed EAM potential was obtained from the literature (Foiles et al., 1986). The normal vibrational analysis was initially performed using a cubic Brillouine zone (−2π/a_{0}, 2π/a_{0}) (a_{0} is the f.c.c. lattice constant) by sampling about 10^{6} phonons. This yielded and in V_{eff}. The classicallike NPT (constant pressure and temperature in a closed system) Monte Carlo (MC) calculations were subsequently preformed using V_{eff} instead of the bare potential for 256 Cu atoms (4^{3} f.c.c. unit cells). Threedimensional periodicity was imposed. 20000 MC steps were calculated to reach the equilibrium and a further 10000 MC steps were performed to obtain the cumulants; each MC step contained 256 movements of atoms and one variation of the lattice constant.
For the thin Cu film, there exists no threedimensional periodicity, thus leading to difficulties in applying the above formula. In a practical sense, however, for the estimation of the V_{eff} can be replaced by the bulk values, although rather crudely from a theoretical point of view. This would give some overestimation of the vibrational amplitude of surface phonons at low temperature since the eigenfrequencies should be smaller than the bulk ones. Within this approximation, the PIECP calculations can similarly be performed. A sixlayer Cu(111) film was assumed, each layer containing 48 atoms in a rectangular lattice; the lowest layer was assumed to be vibrationally fixed. Similar NPT MC simulations were performed to obtain the cumulants of intra and interlayer atom pairs.
cumulants for comparison with the experimental data, the quantummechanical correction in4. Results and discussion
The calculated , together with the experimental observations. Although the calculated values of the distance R are shifted by ∼0.01 Å compared with the and Xray diffraction experimental data, the temperature dependence of R (slope of the curves) indicates that the thermal expansions are comparable. The deviations of the absolute values are partly derived from the fact that the EAM parameters were determined in a classical scheme. Other temperaturedependent quantities, such as the meansquare relative displacement C_{2} and the meancubic relative displacement C_{3}, are also in good agreement with the experimental values. In the C_{2} plot, the Debye fit was evaluated by using the correlated Debye model (Beni & Platzman, 1976) with the Debye temperature as a fitting parameter. It can be concluded here that the PIECP calculations combined with the EAM is sufficiently adequate to describe the thermal properties of bulk Cu.
cumulants for the firstnearestneighbour Cu–Cu shell of bulk Cu are shown in Fig. 1Let us now turn to the results for the Cu(111) film. The temperature dependences of R, C_{2} and C_{3} for the firstnearestneighbour Cu–Cu shells are depicted in Figs. 2, 3 and 4, respectively. In Figs. 2, 3 and 4, the obtained quantities for each layer are given separately. For instance, the curves denoted 12 correspond to the firstandsecond interlayer quantities, while those denoted 11 correspond to the firstandfirst intralayer quantities. Note that the firstandsecond layer (12) distance in Fig. 2, as well as all the intralayer distances 11, 22, 33 and 44, are contracted compared with the bulk one, while the secondandthird (23) and thirdandfourth (34) layer distances are elongated. (Fig. 2), C_{2} (Fig. 3) and C_{3} (Fig. 4) are found to be significantly enhanced for the firstandsecond layer.
We will next compare the present theoretical results with experimental data (Kiguchi et al., 2000). In the previous work, polarizationdependent Cu (and also Ni) Kedge spectra were recorded for ultrathin Cu(111) [or Ni(111)] films grown epitaxically on HOPG (highly oriented pyrolitic graphite) by varying the thickness. By using the results for the thickness dependence of anisotropic C_{2} and C_{3}, the surface outofplane and inplane components can be extracted. Moreover, by assuming (roughly) that all the contributions, other than the inplane firstandfirst (11) and the outofplane firstandsecond (12) contributions, are equal to the bulk contribution, three components of 11, 12 and bulk were successfully determined. For detailed procedures of the data analysis, see the work of Kiguchi et al. (2000). The results are summarized in Table 1. In the experimental work, the parameters were obtained as differences between 100 and 300 K. For C_{2}, the ΔC_{2} values are replaced by the effective Debye temperature θ_{D}, which is similarly evaluated by using the correlated Debye model (Beni & Platzman, 1976).

Although the difference between the outofplane and the inplane vibrations is underestimated in the calculations, the calculated and experimental results agree well semiquantitatively. It is concluded, both experimentally and theoretically, that the surface outofplane vibration is significantly softer and more anharmonic than the inner, lateral or bulk vibration. As a reason for the underestimation in the calculations, one can suppose that the present theoretical calculations are valid for perfect films, while the actual films should contain many defects and some roughening might already have occurred at room temperature. This is indicated by the experimental findings of smaller C_{2} (Kiguchi et al., 2000). Since the surface area becomes wider in the presence of defects, the surface Debye temperature would effectively be lowered. One should also note that the surface Debye temperature of Cu(100) determined by LEED (low energy electron diffraction) is 235 K (Müller et al., 1995), which is still lower than the result of 262 K.
and largerIt is essential to recall intrinsic differences between C_{2} and C_{3} are obtained uniquely by It is natural that the outofplane vibration is enhanced at the surface as long as the absolute displacement is discussed. This is ascribed simply to a lack of atoms above the surface. On the contrary, from a local point of view, the outofplane motion is not always enhanced, as can be seen in the literature: in the p4g(2 × 2)N/Ni(001) surface, the outofplane N–Ni bond has been found to be much stiffer and less anharmonic than the lateral bonds (Wenzel et al., 1990); similarly the outofplane S–Ni bond in c(2 × 2)S/Ni(110) has been found to be stiffer and less anharmonic than the lateral bonds (Yokoyama et al., 1994). In the present case, the surface outofplane bond is found to be weaker and more anharmonic than the inplane bond.
and diffraction. Diffraction techniques such as LEED provide information on absolute displacements with respect to the lattice, while gives relative displacements between Xrayabsorbing and neighbouring atoms. Although within the correlated Debye model the value of the Debye temperature should be equal to that given by diffraction,The aim of this work is to obtain a hint of surface melting, which is an initial stage of bulk melting. The enhancements of the vibrational amplitude and anharmonicity along the surface normal give rise to an important conclusion, since a roughening transition should occur through the `hopping' of surface atoms, for which not only the amplitude but the anharmonicity of the vibration are essential.
5. Conclusions
In this article, the usefulness of the PIECP theory for the evaluation of the
cumulants or, in general, the radial distribution function is demonstrated. In the PIECP theory, the physical quantities converge to be harmonic at low temperature and to be classical at high temperature. Therefore, anharmonicity at low temperature cannot be estimated. In order to obtain such information, one has to perform quantumstatistical perturbation calculations or more sophisticated pathintegral Monte Carlo simulations. Lowtemperature anharmonicity is usually not very important for researchers and the PIECP theory is thus practically very useful to simulate the cumulants. In summary, the PIECP theory can handle the quantum effect, anharmonicity, many degrees of vibrational freedom, threedimensional periodicity, highernearestneighbour interactions and manybody interactions. It is essentially difficult or impossible to include all these contributions using other theories, at least in a practical sense of numerical calculations. It is worth noting again that the EAM matches the PIECP theory excellently, this allowing one to investigate the vibrational properties of metals quantum statistically.In the simulations of thin Cu(111) film, thermal vibration and local
are larger than in the bulk metal. The relative motions, focussing on the local bonds of the surface, are found to be enhanced in the surface normal direction. These findings are consistent with the experiments. The enhanced vibrational amplitude and anharmonicity in the surface normal direction could be the trigger of a roughening transition, surface melting and consequent bulk melting.Acknowledgements
The author is grateful for the financial support of the GrantinAid for Scientific Research (No. 09640598).
References
Beni, G. & Platzman, P. M. (1976). Phys. Rev. B, 14, 1514–1518. CrossRef CAS Web of Science
Cuccoli, A., Giachetti, R., Tognetti, V., Vaia, R. & Verrucchi, P. (1995). J. Phys. Condens. Matter, 7, 7891–7938. CrossRef CAS Web of Science
Cuccoli, A., Macchi, M., Neumann, M., Tognetti, V. & Vaia, R. (1992). Phys. Rev. B, 45, 2088–2096. CrossRef Web of Science
Foiles, S. M. (1985). Phys. Rev. B, 32, 3409–3415, 7685–7693.
Foiles, S. M., Basks, M. I. & Daw, M. S. (1986). Phys. Rev. B, 33, 7983–7991. CrossRef CAS Web of Science
Fujikawa, T. & Miyanaga, T. (1993). J. Phys. Soc. Jpn, 62, 4108–4122. CrossRef CAS Web of Science
Fujikawa, T., Miyanaga, T. & Suzuki, T. (1997). J. Phys. Soc. Jpn, 66, 2897–2906. CrossRef CAS Web of Science
Kiguchi, M., Yokoyama, T., Matsumura, D., Endo, O., Kondoh, H. & Ohta, T. (2000). Phys. Rev. B, 61, 14020–14027. Web of Science CrossRef CAS
Kleinert, H. (1995). Path Integrals in Quantum Mechanics, Statistics and Polymer Physics. Singapore: World Scientific.
Müller, S., Kinne, A., Kottcke, M., Metzler, R., Bayer, P., Hammer, L. & Heinz, K. (1995). Phys. Rev. Lett. 75, 2859–2862. CrossRef PubMed Web of Science
Wenzel, L., Arvanitis, D., Rabus, H., Lederer, T., Baberschke, K. & Comelli, G. (1990). Phys. Rev. Lett. 64, 1765–1768. CrossRef PubMed CAS Web of Science
Yokoyama, T. (1998). Phys. Rev. B, 57, 3423–3432. Web of Science CrossRef CAS
Yokoyama, T. (1999). J. Synchrotron Rad. 6, 323–325. Web of Science CrossRef CAS IUCr Journals
Yokoyama, T., Hamamatsu, H., Kitajima, Y., Takata, Y., Yagi, S. & Ohta, T. (1994). Surf. Sci. 313, 197–208. CrossRef CAS Web of Science
Yokoyama, T., Ohta, T. & Sato, H. (1997). Phys. Rev. B, 55, 11320–11329. CrossRef CAS Web of Science
Yokoyama, T., Yonamoto, Y., Ohta, T. & Ugawa, A. (1996). Phys. Rev. B, 54, 6921–6928. CrossRef CAS Web of Science
© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.