Projecting picosecond lattice dynamics through x-ray topography

Projecting picosecond lattice dynamics through x-ray topography

A method for time-resolved x-ray diffraction studies has been demonstrated. As a test case, coherent acoustic phonon propagation into crystalline InSb is observed using a laser plasma x-ray source. An extended x-ray topogram of the semiconductor’s surface was projected onto a high spatial resolution x-ray detector and acoustic phonons were excited by rapidly heating the crystal’s surface with a femtosecond laser pulse. A correlation between the spatial position on the x-ray detector and the time of arrival of the laser pulse was encoded into the experimental geometry by tilting the incident laser pulse with an optical grating. This approach enabled a temporal window of 200 ps to be sampled in a single topogram, thereby negating the disadvantages of pulse-to-pulse fluctuations in the intensity and spectrum of the laser-plasma source.

Both laser plasma x-ray sources and synchrotron radiation facilitate time resolved x-ray diffraction studies of light induced structural changes on a time scale of picoseconds. Notable successes at synchrotron sources have followed light induced structural rearrangements in crystals of macromolecules,1,2 have enabled the observation of coherent acoustic phonons,3,4 and have quantified large amplitude motions within a photoexcited chemical system.5 Laser plasma x-ray sources have facilitated gas phase6 and solution phase7 picosecond x-ray absorption experiments, and have enabled the growth of disorder in an organic sample,8 the nonthermal melting of a semiconductor surface,9 and the consequent propagation of a strainwave through a semiconductor lattice,10 to be visualized at near ~or even sub! picosecond
resolution.

When visualizing rapid structural rearrangements, the very high peak x-ray brilliance of synchrotron radiation provides a tremendous advantage. The choice either a polychromatic or a tunable monochromatic x-ray source with low angular divergence ~as well as low pulse-to-pulse fluctuations! enables considerable freedom when designing an experiment. The primary disadvantage of synchrotron radiation is the relatively long duration of the x-ray pulses, which cannot be pushed significantly below 50 ps without severe loss of x-ray intensity.11 The use of cross-correlation techniques12 or x-ray streak cameras,4 however, have enabled improvements in the temporal resolution beyond the pulse length of synchrotron radiation when studying the laser induced melting and reordering of semiconductor surfaces. Schemes for generating synchrotron based x-ray pulses with a duration of a few hundred femtoseconds have also been pursued,13,14 and used for time-resolved x-ray diffraction
studies.15

Laser plasma x-ray sources provide pulses approximately two10 to three8,9 orders of magnitude shorter in duration than those available at synchrotrons. Although the peak number of x-ray photons generated is comparable with that generated at a synchrotron source, experiments using a laserplasma source are hampered by a high beam divergence and large pulse-to-pulse fluctuations. In particular, laser plasma generated x-ray photons are emitted over the full sphere ~108 photons/pulse in this study!, and the integrated intensity can fluctuate by one or two orders of magnitude. Nevertheless, by accumulating a large body of experimental data, Rose Petruck et al.10 succeeded in piecing together a movie following strain-pulse propagation through a semiconductor lattice after the surface was rapidly heated by a femtosecond ~fs! laser pulse. The data were of sufficient quality that a quantitative description of the lattice dynamics could be given over a 250 ps window.

In this work an experimental configuration is presented which simultaneously takes advantage of natural beam divergence and negates problems associated with pulse-to-pulse fluctuations of a laser plasma x-ray source. The experimental setup is illustrated in Fig. 1 and is motivated by a proposal to use topography as a methodology for recovering subpicosecond temporal resolution in x-ray diffraction experiments.16 A crystalline sample of InSb was first aligned relative to the iron target laser plasma x-ray source so as to satisfy Bragg’s condition for diffraction from a surface. Due to the natural divergence of the x-ray beam, an x-ray topogram could be recorded on a high-spatial resolution charge coupled device ~CCD! camera. A 100 fs laser pulse ~l5800 nm, 2.5 mJ/ pulse! was used to rapidly heat the crystal’s surface ~experimental area '25 mm2! and its arrival relative to the x-ray

pulse was set by an optical delay line ~not illustrated in the schematic!. Although a single laser pulse heated the entire surface exposed to x rays, an additional time delay at different spatial locations along the crystal’s surface was introduced by reflecting the laser pulse from a diffraction grating prior to it being imaged onto the sample ~Fig. 1!. As such the reflected wave front became tilted relative to its direction of propagation, with the result that different rays traversed different optical path lengths before arriving at the semiconductor’s surface. Since the maximum difference in the optical path lengths of different rays was 6 cm, a temporal window of 200 ps was sampled in a single x-ray topogram ~Fig. 1!, reflecting the finite time required for the tilted wave front to sweep across the crystal’s surface exposed to x rays.16

Through this approach the disadvantages of pulse-to-pulse fluctuations in the intensity and spectrum of a laser-plasma source were negated, since all sampled time points were exposed to an identical ~accumulated! x-ray probe and, hence, were directly comparable. Figure[ 2a] illustrates a typical x-ray topogram from a crystal of InSb without any laser induced heating, as recorded on a CCD detector with high spatial resolution 25 mm per pixel!. Because the x-ray source is polychromatic, different wavelengths satisfy the Bragg diffraction condition at different angles of incidence, hence, the two Fe Ka lines of the laser plasma’s target are projected as two parallel lines on the x-ray topogram [Fig. 2a]. When the sample is rapidly heated by the fs laser pulse, a structure in the x-ray topogram becomes visible [Fig. 2b], with one Fe Ka line first being split, then remerging, as the position along the detector’s x axis is increased. This change in the x-ray topogram reflects the fact that laser induced heating first perturbs the crystalline lattice and later, as energy disperses throughout the crystal, the original lattice spacing is recovered. The dependence of the measured scattering intensity on the spatial position of the detector is due to the correlation between the spatial position and the time of arrival of the laser pulse, which was encoded into the experimental geometry [Fig. 1].This results in an image which is a time-resolved rocking curve.

laser pulse was altered by changing the optical delay line, all of the laser induced features of the x-ray topogram [Fig. 2b] were reproduced, but were translated along the x axis of the detector. This observation is quantified in Fig. 3, which plots the translation of the x-ray topogram versus the temporal change in the optical delay line. From this result the effective temporal resolution of the x-ray detector could be calibrated as 1.2 ps per pixel, which is in agreement with that determined from the optical configuration.

Although the acquisition time required to sample this 200 ps window was relatively short data shown in [Fig. 2b]derived from a single 60 min run the diffraction data were of sufficiently high quality that a good theoretical fit to the experimental results could be obtained. In [Fig. 2c] we show a simulated x-ray topogram for comparison. A detailed description of the theoretical model is given in Refs. 17 and 18

and it has proven to be applicable to a broad spectrum of experimental results.3,10 Intuitively, energy from the fs laser pump couples into the semiconductor by exciting electrons from the valence to the conduction band. After absorption, energy is initially dispersed as optical phonons, which are characterized by specific motions of atoms within the unit cell. On a slower time scale coherently excited acoustic phonons are generated,3,10 which are characterized by specific motions of the crystalline lattice and their propagation is determined by the optical, electronic and acoustic properties of the material.17 The propagation of these acoustic waves through the semiconductor perturbs the lattice spacing on a time scale of picoseconds, and results in the splitting and subsequent relaxation! of the Fe Ka line modeled in [Fig. 2c].

Several extensions to the experimental protocol presented here are immediately foreseeable. The use of a large area x-ray camera with a high spatial resolution would enable several diffraction spots to be recorded simultaneously, thereby recovering more structural information than that accessible from surface diffraction alone. Furthermore, the potential temporal resolution can be manipulated by the choice of experimental geometry. Should the fs laser pulse sweep directly across the sample’s surface16 then the effective temporal resolution of each pixel of the x-ray camera would be a few tens of femtoseconds. Such a change in the experimental geometry would not result in any loss in signal-to-noise, although the temporal resolution achievable in practice, even using a deconvolution algorithm, cannot be pushed significantly beyond the duration of the x-ray probe.19 Nevertheless, when used in conjunction with sufficiently short x-ray pulses, the experimental approach presented here promises a direct route to observing subpicosecond phase transitions during laser induced melting of well-diffracting samples.

The authors thank Michael Wulff and Friedrich Schotte for assistance with preparatory experiments at ID09 of the European Synchrotron Radiation Facility [ESRF]. The authors acknowledge support from the European commission ‘‘Improving the Human Potential program.’’ J.L. and R.N. acknowledge support from the Swedish Science Research Council [V.R.]and the Swedish Strategic Research Foundation (S.S.F)

Expand+Unusual lattice dynamics of vanadium under high pressure

Unusual lattice dynamics of vanadium under high pressure

The electronic structures and lattice dynamics of pressure-induced complex phase transitions [bcc → hR1(110.5°) → distorted-hR1(108.2°) → bcc] in vanadium as a function of pressure up to 400 GPa have been investigated with an ab initio method using density functional perturbation theory (DFPT). At ambient pressure, the soft transverse acoustic phonon mode corresponding to Kohn anomaly appears at a wave vector q = 2k F along [ξ00] Γ → H high symmetry direction. The nondegenerate transverse acoustic branches TA1 on 〈110〉 and TA2 on 〈001〉 show an exceptionally large split at high symmetry point N (0.5 0.5 0.0). The lattice dynamical instability starts at a pressure of 62 GPa (V/V 0 = 0.78, where V 0 is experimental volume of bcc-V at ambient conditions), derived by phonon softening that results in phase transition of bcc → hR1 (alpha = 110.5°). At compression around 130 GPa (V/V 0 = 0.67), the rhombohedral angle of hR1 phase changed to 108.2°, and the electronic structure changed drastically. At even higher pressure, ≈250 GPa (V/V 0 = 0.57), lattice dynamic calculations show that the bcc structure becomes stable again.

A vast number of extremely interesting and fundamental physical and chemical changes in solids have been reported in high-pressure experiments. The Kohn anomaly in bcc transition metals is considered to be one of the unsolved problems in the field of high-pressure structural phase transition. Kohn anomaly is characterized by electron–phonon features. Vanadium is an interesting case because this metal is one of the so-called hard superconductors (1) and it has been shown that the application of high pressure increases the superconducting temperature Tc (2, 3). The experimental phonon dispersion curves at ambient conditions in bcc phase using inelastic neutron scattering has been measured and for the longitudinal (111) branch it exhibits a dip near the (2/3 2/3 2/3) position (4). In another experiment by Page (5), a peak in the frequency distribution for vanadium is observed at 2 × 1012 c/s, and this peak is attributed to Kohn anomalies (5). As far as calculations are concerned, the lattice dynamics study of bcc vanadium at high pressure using a full-potential linear muffin-tin orbital (LMTO) method has been performed by Suzuki and Otani (3), and these calculations show a strong change of phonon density of states with volume at 130 GPa. It has been found that the transverse acoustic phonon mode TA [ξ00] around ξ = 1/4 shows a dramatic softening under pressure and becomes imaginary at pressures >130 GPa, indicating the possibility of a structural phase transition. Landa et al. (6) have performed the elastic constants calculations for vanadium using exact muffin-tin orbital (EMTO) and shown the trigonal shear elastic constant (C44) softening at ≈200 GPa due to the intraband nesting of the Fermi surface (FS).
Very recently, Ding et al. (7) performed the high-pressure experiment for bcc vanadium using diamond-anvil cell and synchrotron x-ray diffraction and, for the first time, observed a transition from bcc to a rhombohedral phase (hR1) at 69 GPa; this is a new type of high-pressure structural phase transition that was not found in earlier experiments for any element or compound. Therefore, it will be very interesting to perform first-principles calculation using this new high-pressure phase to look for the origin of such structural phase transitions. Here, we present calculations for the phonon dispersion, phonon density of states, and electronic band structures for different phases as a function of pressure. Details about the performed calculations can be found in Methods. We observed the structural sequence bcc → hR1(110.5°) → distortion-hR1(108.2°) → bcc with increasing pressure.

Results and Discussion

We performed the calculations for the phonon dispersion curves along the high symmetry direction Γ–H–P–Γ–N for bcc-V at ambient conditions (Fig. 1). The solid lines are obtained by calculation, and the open symbols are the experimental data taken from the experimental work that were observed with x-ray thermal scattering (4). There is a steep slope observed at 1/3(Γ–H), caused by Kohn anomaly, on computed transverse acoustic mode (TA) along [ξ00] direction. Our calculations are in overall good agreement with experimental results, including anomalous frequency distribution near high-symmetry point N along Γ–N. There are two branches of one LA mode and other degenerated TA1 (〈1̄10〉) and TA2 (〈001〉) modes along the [ξξ0] direction. There exists a crossover of LA mode with TA2 mode at zone boundary point N. The transverse dispersion curve for [ξξ0] TA2 is higher than those of the longitudinal dispersion curve near the BZ of [ξξ0] symmetry direction. The present computed result gives better ordering of the branches than previous calculations (8, 9).

Fig. 2 shows the phonon density of state (PDOS) for bcc-vanadium (bcc-V) as a function of pressure. The lower frequency phonon mode caused by the Kohn anomaly is found to be grown and softened, but all other phonon frequencies move to higher values with increasing pressure. At V/V 0 ≈ 0.78 (62 GPa), negative frequency appears in the PDOS spectrum, but the PDOS still shows the typical bcc structural spectrum. After V/V 0 = 0.78 (62 GPa), the PDOS shows a

different profile of frequency distribution, indicating the dynamic instability of bcc-V at V/V 0 ≈ 0.78.

After comparing our calculated phonon spectrum with experiments performed under ambient conditions, we examined the effect of pressure on the phonon spectrum. Here, we chose three volumes: V/V 0 = 0.78, 0.75, and 0.56. Fig. 3 a shows that, at values up to V/V 0 = 0.78, bcc phase with the phonon mode linked to a Kohn anomaly has softened to negative values, and this can induce a high-symmetry loss and breaking in the H direction. This also drives the distortion in the bcc phase along (100) direction. The phonon spectra near Γ point also show softening. At V/V 0 ≈ 0.75 (112 GPa), the hR1 structure with α = 110.5° is more stable in our calculations and the corresponding phonon dispersion at this pressure is presented along high-symmetry direction Γ–FA–T–Γ–L (Fig. 3 b). As pressure increased to 130 GPa, i.e., volume compression V/V 0 = 0.67, the rhombohedral angle of hR1 phase changes to 108.2°. At V/V 0 = 0.56 (260 GPa), the rhombohedral angle changes back to 109.47° (i.e., the bcc phase), indicating that the bcc phase is dynamic stable or hR1 with an ideal angle of 109.47°. The phonon dispersion at this pressure is calculated and shown in Fig. 3 c. Compared to the ambient condition bcc phonon dispersion spectrum, there appears a deep dip of LA mode in (ξξξ) direction. This phonon softening in high-pressure bcc phase might make bcc phase unstable under extreme conditions (e.g., >500 GPa). Although the phonon dispersion of high-pressure bcc phase indicates dynamic stability at this pressure, the Kohn anomaly is still there. The anomaly in frequency at symmetry point N is absent in the high-pressure bcc phase.

The hR1 structure of V can be easily characterized starting from a parent bcc structure with primitive unit cell and considering the atomic displacements along the cubic diagonal, which transform as a transverse zone boundary instability at point N of the bcc Brillouin zone (BZ). At ambient pressure, vanadium crystallizes in a body-centered cubic phase Im3̄m, which has anomaly at symmetry point N mainly due to the larger elastic anisotropy (10). This is why, in our studies, the frequencies of transverse branch at N point are also found to be higher than those of longitudinal branch. We have shown that the frequencies of longitudinal mode LA and the transverse mode TA2 [polarized along (001)] mode as a function of pressure (Fig. 4). The frequency dispersion of LA and TA2 at V/V 0 ≈ 0.80 change back to the normal order (i.e., TA2 mode frequency is lower than LA mode frequency).

To understand high-pressure structural phase transition of vanadium by dynamical stability, we plotted the frequencies at Γ center along high-symmetry Γ–H and Γ–N directions of bcc BZ (Fig. 5, solid lines). The frequencies at Kohn anomaly [H(1/3 0 0)] soften mode are shown by dashed lines as a function of pressure. As mentioned above, the soft modes caused by Kohn anomalies are in the lower-frequency part along Γ–H direction. Under compression, the frequencies at (1/3 0 0) soften, and the phonon modes smoothly and slowly go to negative values up to V/V 0 = 0.79. It is reasonable to suggest that the phonon mode softening related to Kohn anomaly drives the structural phase transition from bcc to hR1. Above this pressure (up to V/V 0 = 0.59), the phonon dispersion of bcc phase shows a complete dynamical instability. At Γ point, all of the phonon branches along high-symmetry directions are negative, suggesting that no active phonon modes are present in bcc phase. This finding also indicates that bcc symmetry is broken under high pressure. It is noteworthy that the angle of rhombohedral hR1 phase changes from 110.5° to 108.2° at V/V 0 ≈ 0.67. A further compression changes this angle back to 109.47° at V/V 0 = 0.57. The origin of these two changes in rhombohedral angle α cannot be explained by the lattice dynamic calculation alone. Unlike the phase transition of bcc to hR1 at V/V 0 = 0.78, the transition of hR1 with angle 108.2° to bcc at V/V 0 = 0.57 shows an abrupt change in phonon frequency.

At ambient pressure, vanadium crystallizes in a body-centered cubic phase Im3m, and bcc structure can also be shown in rhombohedral structure. If the rhombohedral angle is 109.47°, then this structure will be equal to the bcc structure. To avoid the effect of symmetry, we have done all of the electron band structural calculations for vanadium using rhombohedral structural setting. The electronic band structure of bcc and hR1 phases as function of volume (or pressure) are shown in Fig. 6 a–d. For comparison, we have also plotted the band structure of bcc (Im3m) phase at ambient condition (Fig. 6 e).

From ambient volume to V/V 0 = 0.80, no pronounced change can be seen from the electronic band structure. In Fig. 6 a, we show the electronic band structure at volume compressed to V/V 0 = 0.80. There are two conduction bands: band 1 has symbols Γ1–FA1–T1–L1 and band 2 has the symbols Γ2–FA2–T2–L2 at high symmetry points, respectively. Both bands are partly filled by conduction electrons and above Fermi level at Γ point. In the rhombohedral hR1 structure with a 110.5° angle, shown in Fig. 6 b, there are still two conduction bands. The noticeable difference in hR1 band structure compared to band structure in Fig. 6 a is that, at Γ point, band 1 moves below the Fermi level and separates from the other two bands. Under pressure (i.e., at V/V 0 = 0.70), the Fermi level moves up and conduction band 1 is almost full filled by conduction electrons (data not shown). There is only one band that crosses the Fermi level, just like the electronic band structure of fcc metals.

At V/V 0 = 0.65, the rhombohedral angle α changes from 110.5° to 108.2°; the corresponding band structure is shown in Fig. 6 c. For α = 108.2° rhombohedral phase, band 2 moves from above to below the Fermi level and links with band 1 at Γ point, as shown in Fig. 6 a. At this point, rhombohedral angle α transforms from 110.5° to 108.2°; it is reasonable to suggest that this transition is driven by the electronic phase transition shown in Fig. 6 c. At even higher pressure, in high-pressure bcc phase, two conduction bands connected below the Fermi level at Γ center as shown in Fig. 6 d. One can see a huge change in electronic band structure going from Fig. 6 c to d. An empty band (≈1 eV above the Fermi level at point Γ) at V/V 0 = 0.60, as shown in Fig. 6 c, moves down as pressure increases and merges with two occupied bands at point Γ (Fig. 6 d). This change in electronic structure is responsible for an abrupt change in phonon frequencies at the transition from hR1 (108.2°) to high-pressure bcc phase.

We have also calculated the bulk modulus for vanadium using the Birch–Murnaghan equation of state for ambient bcc and hR1 phases. Our calculation gives the bulk modulus for bcc phase, B = 177.7 GPa with B′ = 4.3, and hR1 phase, B = 182.5 GPa with B′ = 3.6. These data compare very well with our experimental results (7).

In conclusion, we have performed ab initio calculations of lattice dynamics and electronic band structure as a function of pressure for vanadium. Based on our lattice dynamics calculations, we show the following structural phase transitions in vanadium bcc → hR1(110.5°) → distorted-hR1(108.2°) → bcc as a function of pressure. Our phonon calculations confirm the phase transition from BCC to hR1, which takes place at ≈62 GPa and compares very well with our experimental value of 69 GPa (7). At high pressure, based on our phonon calculations, we predict two new phase transitions: one from hR1 (110.5°) to distorted-hR1 (108.2°) at ≈120 GPa, and finally the reentrance of bcc phase at ≈250 GPa. We have explained our lattice dynamics results in terms of electronic band structure. A huge change in the electronic structure is driving force behind the structural phase transition. We welcome more experiments to look into the unusual lattice dynamics of vanadium.

Methods

Electronic structural properties and phonon frequencies in vanadium are analyzed by using the density-functional perturbation theory (DFPT) (11). The vibrational properties have been computed by using a recent implementation of ultrasoft pseudopotentials (12) into DFPT formalism in the PWSCF (Plane-Wave Self-Consistent Field) code (www.pwscf.org). We have used Vanderbilt ultrasoft pseudopotential with 3s, 3p, 3d, and 4s states as valence states (13). The calculations are carried out by using generalized gradient approximation for the exchange and correlation potential as described by Perdew et al. (14).

The electronic wave functions are represented in a plane-wave basis set with a kinetic energy cutoff of 50 Ry. The Brillouin zone (BZ) integrations are carried out by the Gaussian smearing technique using a 12 × 12 × 12 k point mesh with shift from origin. Total energies are converged to within 10−4 Ry/atom with respect to energy cutoff changed from 30 to 90 Ry. The dynamical matrices were calculated on a 6 × 6 × 6 grid. All of the structural parameters are fully relaxed

COMPARISON OF THE LATTICE-DYNAMICS

COMPARISON OF THE LATTICE-DYNAMICS
THERMODYNAMIC PROPERTIES

Several paths are available for calculating macroscopic thermodynamic properties for a system with given interparticle forces. Most paths are approximate. The lattice dynamics approximationl) treats correctly all terms in the energy which are quadratic in the particle displacements. The cell-model approximation2) in which a single particle moves in the field of its fixed neighbors modifies the quadratic terms but includes an estimate of anharmonic' corrections. More sophisticated theories3) treat anharmonic perturbations analytically (the actual calculations require fast computers) but seem complicated enough to attract few follo\vers.

The approximate methods have the important advantage of being quick and inexpensive to calculate. The computer experiments giving exact thermodynamic properties, either by following the motion of the particles or by sampling the configuration space 4), are relatively expensive because they require so much computer time, particularly if high precision is necessary. To halve the statistical uncertainty in computer-experiment results requires quadrupling the computer time used up. The precision also depends on the sensitivity of the property measured to fluctuations in pressure and energy. Because successively higher derivatives of the free energy involve higher moments of the pressure-tensor component distributions and energy distributions, the time necessary to characterize derivatives increases rapidly with derivative order. Second-order elastic constants and the specific heat involve second moments; third-order elastic constants involve third moments; and so on. So far only first-and second-derivative quantities have been examined. Thus, if the approximations should prove to come close to exact results, they would give us a useful shortcut to accurate thermodynamic properties. It was in the hope of establishing their usefulness that we undertook these calculations.

In this paper we compare the results of solid-phase Monte Carlo experiments on 108 particles interacting with the Lennard-Jones and exponential-six potentials with the predictions of 108-particle lattice dynamics and the cell model. In addition to the energy and pressure, we compare the secondderivative quantities: specific heat, Griineisen y, and elastic constants, with approximate predictions. Although either computer method, molecular dynamics or Monte Carlo, can be extended to quantum systems by using the \Vigner-Kirkwood Planck's constant expansion of the free energy5), we have made our comparisons using classical mechanics. To find out how important the small size and classical nature of our systems are, we use the latticedynamics method to investigate the number dependence of all of the thermodynamic properties and quantum corrections to the elastic constants. Quan-t tum corrections to other thermodynamic properties have been calculated elsewhere6).

In section 2 we describe the lattice-dynamics calculations. The method has been in use for over 50 years, although many so-called lattice-dynamics calculations of elastic constants have actually been calculations of the elastic response of a static lattice. We take lattice vibrations explicitly into account. The first correct harmonic calculation using normal-mode vibrations was announced by Feldman7). He obtained expressions for the second-order elastic constants which involved the vibration frequencies and their strain derivatives. The frequencies were then obtained by the usual method of diagonalizing the dynamical matrix. The frequency-shift derivatives were calculated by means of perturbation theory. We calculate the work of deforming the crystal by the alternative procedure of computing the free energy numerically for several slightly different values of the strain and then fitting the results to a strain polynomial. Besides avoiding the tedious algebra of Feldman's analytic approach, the numerical method is more readily generalized to higher-order elastic constants.

In section 3 we describe the cell-model calculations. The cell model, although actually only appropriate for solids, was first used in an attempt to describe gases and liquids 8). In the cell model the effect of heating the crystal lattice is approximated by a "one-particle" model in which a single particle moves in the field of its fixed neighbors. Neglecting interparticle correlations by approximating an N-body problem by a one-body problem for classical systems, most reasonable at low temperatures and is exact only in the static-lattice limit.

The advantage of the cell model over the harmonic approximation lies in its ability to estimate anharmonic contributions from potential-energy terms beyond the quadratic ones. The cell model has often been used to calculate energy, pressure, and specific heat 9). Our elastic-constant calculations are a new use of this model. Henkel10) has studied a similar model, a quantum cell model in which the potential was expanded in powers of displacement and the quartic contributions were treated by perturbation theory. If the perturbations were ignored, Henkel's work would reduce to the usual harmonic Einstein model.

In section 4 we compare the tabulated results from both approximations and consider the dependence of the results on number of particles. Vie also discuss some interesting cancellations found in the course of the Monte Carlo calculations. In section 5 we assess the importance of quantum effects on the elastic constants.

2. Lattice-dynamics calculations.

The lattice dynamics calculations are based on the approximation of truncating a Taylor expansion of the lattice potential energy after the quadratic terms in the particle displacements. Sometimes this is called the "quasi-harmonic" approximation. The coefficients in the Taylor series expansion are calculated from the assumed force law, and the lattice sites are chosen to match the structure of the lattice being described. The expansion is made about a configuration in which each atom is fixed at its average position in a perfect crystal with fixed center of mass. Thus the truncated potential depends upon the size and shape of the assumed static lattice configuration. By changing to normal-mode coordinates the quasi-harmonic Hamiltonian can be rewritten as a sum of 3N 3 independent harmonic oscillator Hamiltonians. The partition function of an oscillator, either quantum or classical, is known11) so that the quasi-harmonic thermodynamic properties of the system can be calculated.

The partition function, Z = exp(-A/kT), where A is the Helmholtz energy and kT is Boltzmann's constant times the absolute temperature, can be written as a product of single oscillator partition functions:

In the classical limit only the first term in the expansion, kT/hv, is kept. from the quadratic Hamiltonian; the center of mass contribution is Zcm. For accurate work on small crystals Zcm has to be included if comparisons are made with Monte Carlo calculations in which the center of mass is allowed to move.

Thermodynamic properties can all be derived from the partition function. Temperature derivatives can be evaluated explicitly to compute the energy and specific heat:

"Strain" derivatives are harder to evaluate. The strains are defined in terms of three vectors colinear with the edges common to one corner of a parallelepiped produced by deforming a cube of crystal with sidelength a. If aI, a2, and a3 are the vectors, then

are the six independent strains. The dynamical matrix then gives a complicated implicit relation for the vibration frequencies as functions of the strains. Because there is no convenient expression for v(1]) it is easiest to proceed numerically. Both the average pressure tensor component

where fi and ri are the Cartesian coordinates of rand r, and 1]1 is the only nonzero strain. In the strained configuration the Hamiltonian is again expanded, the quadratic terms kept, and the result diagonalized, giving a new set of frequencies and the free energy A (1]1). The diagonalization can be visualized in terms of plane-wave solutions of the classical equations of motion. The waves propagate through the crystal with wavelengths and directions imposed by the shape of the crystal and described by wave vectors chosen from a convenient Brillouin zone12). If y is a wave vector in the unstrained crystal then the corresponding y in the strained crystal is:

where Yi and Yi are the Cartesian coordinates of y and y. Applying the Born-von Karman normal-mode analysis for several values of 1]1, the first four or five coefficients in the expahsion,

can be determined. By choosing values of 1J1 separated by 0.0001, both the pressure and Cil were determined with four-figure accuracy in this way. By considering two simultaneous strains along the one axis and the two axis, and analyzing the resulting free energy changes as a series in 1]1 and

The generalization of this technique to calculation of higher-order elastic constants or to mixed strain-temperature derivatives is straightforward. For crystals of lower symmetry one needs to calculate more elastic constants and hence one considers more different combinations of strain. Other strains would also be needed for cubic crystals to determine the six nonzero thirdorder constants or the eleven nonzero fourth-order constants and the mixed strain-temperature derivatives.
The adiabatic elastic constants

where 5 is the entropy of the crystal. Because 5(1]1) and (oT/o5)n can be calculated exactly for a quasi-harmonic crystal, the correction term can be evaluated numerically.
In order to make a comparison with the 108-particle Monte Carlo calculations13, 14) we have calculated the thermodynamic properties for a 108particle system with the same periodic boundaries and Hamiltonian as those used in the Monte Carlo work. Using the nearest-image convention, each particle in the crystal interacts with 107 neighbors according to the LennardJones 6-12 potential

These potentials have both been used principally to describe rare gases and have shown themselves to fit these rea(materials well. Although we picked these potentials because of their value in describing rare gases we expect

that our general conclusions in comparing approximate calculations with exact computer experiments will be valid for potentials describing interactions in salts or metals as well.
Table I gives both the Monte Carlo and the lattice-dynamics results. The three temperatures span the range from about 0.48TtriPIO to O.95TtriPle and the densities correspond closely to zero pressure. Because the static-lattice contributions to the thermodynamic properties present no theoretical problems (they are correctly calculated by any theory) we have tabulated separately the thermal contributions to the thermodynamic properties. The data in the table show that the lattice-dynamics elastic constants are quite· close to the Monte Carlo values at three different temperatures and for both potentials tested.

Structure and lattice dynamics of the ordered phase of solid C70

Structure and lattice dynamics of the ordered phase of solid C70

Since the discovery of a method to produce fullerenes in large amounts [ 1 ] many studies have been performed on these exciting all-carbon molecules. First experimental results on the structure of high purity solid CT0 were presented by Vaughan et al. [ 21 over a year ago. Recently, experimental data on the crystal structure and dynamics of Cl0 single crystals have become available [ 3,4]. The CT0 molecule can be described as two half Cbo molecules joined together by a belt of hexagons [ 5 1, as experimentally confiied by Johnson et al. [ 61 using ZD-NMR. At high temperatures, above 350 K, the elongated-soccerball shaped CT,-m, olecules can rotate quasi freely in the crystal lattice. At these temperatures there is a coexistence of two crystal structures, face centered cubic (fee) and hexagonally close packed (hcp), of which the former one seems to be favoured at still higher temperatures [2]. Crystals grown by sublimation
at 870 K have predominantly the fee growth habit, although a minor fraction of the crystals grown has the hcp morphology. X-ray diffraction [ 3 ] shows that the hcp-grown crystals also exhibit microscopically an ideal hcp structure at high temperatures with rotations around both the long and short molecular axes. Upon cooling these crystals to room temperature the rotations around the short axes freeze out and the long axes of the molecules align along the hexagonal c axis, resulting in a deformed hcp structure.

Somewhat below room temperature the orientational disorder around the long axis freezes out
as well. For this phase an ordered superstructure with four molecules per unit cell is observed. It should be noted that no well defined low temperature crystal structures are obtained when fee-grown crystals are cooled down. In this Letter, harmonic lattice dynamics calculations on the low temperature monoclinic structure of solid C,,, are presented. For the calculations we used rigid CT0 molecules with only van der Waals interaction between them. The structure resulting from a lattice energy minimization is compared to the experimentally determined structure and insight is gained in the orientational ordering of the molecules. The frequencies of the lattice modes are calculated for a variety of interaction potentials.

Lattice dynamics calculations
To perform lattice dynamics calculations on the low temperature monoclinic phase of solid CT0 we applied a standard harmonic lattice dynamics program for molecular crystals, details of which can be found elsewhere [ 71. Rigid CT0 molecules are used, with the atomic positions within each molecule deduced from the five, experimentally determined, nonequivalent atomic positions [ 81 by applying the symmetry operations of the I&s point group. The calculations are performed for zero temperature where an orientational ordering exists, i.e. the molecules do not rotate but have only small amplitude librations around a well-defined equilibrium structure. The interaction between the molecules is modeled as the sum over atom-atom interactions. To describe the atom-atom interaction several potentials known from literature, both of the Lcnnard-Jones and of the expd type, are used #‘. It should be explicitly noted that no Coulomb interaction terms have been included, as will be discussed later. The lattice sum is performed over all molecules within a 25 A range. The lattice energy is minimized with respect to the lattice parameters and with respect to the positional coordinates and the Euler angles of the four molecules in the unit cell. The experimental values are used as the initial values of the 30 minimization parameters. No symmetry restrictions have been imposed on any of these parameters. For the energetically most favourable lattice structure thus obtained the frequencies and eigenvectors of the librational and translational inter-molecular modes are computed.

The results calculated for wavevector R = 0 can be directly compared with experimental Raman data. The transformation properties of the eigenvectors under the symmetry operations of the P2,/m space group enable us to make an unambiguous symmetry assignment of the modes. The experimentally observed structure that is used as starting input for the structure optimization is a monoclinic structure with the unique b,,, axis along the former hexagonal c axis and an angle p close to, but not necessarily identical to, 120”. The a, axis is identical to the hexagonal a axis. It is observed experimentally that the monoclinic unit cell is doubled along the c, direction, which corresponds to one of the former a axes. This doubling of the unit cell is ascribed to an orientational ordering in the former hexagonal plane. It is assumed that in the low temperature structure all CT0 molecules are aligned with their long axes parallel to the b, axis. Given this assumption, which is shown to be valid in the calculations, the orientation of the molecules can be described by the angle @ by which the molecules are rotated around their long axis. At @=O’ the comer points of the top (and bottom) pentagons of the CT0 molecules are directed along the minus h axis. From the experimental data a model describing the orientational ordering was proposed [ 3) in which the molecules in a row along the am axis are oriented identically and molecules in two adjacent rows differ by 180” [ mod.72” ] in orientation. It was argued that this arrangement of molecules would yield the deepest lattice energy minimum.

Results and discussion
The crystal structure obtained via minimization of the lattice energy is in good agreement with the experimentally determined structure. In table 1 the results obtained for various potentials are listed, together with the experimental results. The values for the lattice parameters u, b and c agree to within two percent with the experimentally determined values. The angles LY and y are always found to be 90”. The angle /3 is found to be close to 120”, as can be seen in table 1. These results are in agreement with similar calculations by Guo et al. [ 1 1 ] . The positions of the molecules within the unit cell deviate only slightly ( G 1%) from the ideal lattice positions. The orientations of the long axes of the molecules are found to be parallel to the b, axis, so the calculations justify the assumption mentioned earlier. The energetically most favourable lattice structure is a monoclinic structure, with a doubling of the unit cell along the c, axis due to the orientational ordering, which involves a rotation of the molecules around their long axis. It is important to note that the calculated lattice structure obeys the experimentally observed monoclinic symmetry; the angles #i of molecules i (is 1, . . . . 4) are related via & (~~)=180°[mod.72”]+@2 ($1). The 2,/m screw axis is parallel to the b, axis and intersects the a,+,, plane at position (x, z) with x=0.1663 and y=O.3390, still reflecting the apparent hexagonal symmetry of the high temperature phase. The spacegroup is therefore P2,/m. The molecules are located at positions with site symmetry m. The mirror planes perpendicular to the screw axis

cut through the middle of the & molecules at (0, angle &,, as indicated in fig. 1. This angle @a, as ob- (0, 0) and (0.322, 0.178, 0.500). tained for various potentials, is also given in table 1, As seen from table 1 the calculated lattice energies for the various potentials differ considerably. The aim of the calculations was primarily to understand the CT0 crystal structure and the inter-molecular mode spectrum. These quantities are mainly sensitive to the shape of the interaction potential around the equilibrium position and not that much to the absolute value of the binding energy. No experimental data on the absolute value of the lattice energy for low temperature monoclinic C& are available. Comparing the lattice energy with the experimentally determined high temperature heat of sublimation of 45 kcai/mol [ 12 ] shows that the LJ( 1) potential gives the most realistic lattice energy value. From the calculations some additional conclusions on the orientational ordering can be drawn. The ordering found by the calculation, shown schematically in fig. 1, differs somewhat from the simple model proposed by Verheijen et al. [3] mentioned above. Considering only the P2Jm symmetry two independent angles, & and q&, are needed to unambiguously describe the structure. In the model of Verheijen et al. it is assumed that these angles are $i=O”[mod.72”] and d~=180”[mod.72”]. From the calculations it follows, however, that the difference @i - g2 is not equal to 180” [mod.72” 1. Instead, it is found that the (in principle) independent angles $, and & are related via &=180° [mod.72”] -$,.

Therefore the structure can be described by a single To gain more insight in the nature of the orientational ordering fig. 1 shows a view from the center of molecule 1, in the origin of the unit cell, towards molecule 2. It can be seen that the relative orientation of the molecules is such that the C atoms tend to face the hexagons of the neighbouring molecule. To correctly describe the orientational ordering in solid CbOit is essential that a Coulomb interaction is included in addition to the C-C van der Waals interaction [ 13 1, In a Cbo molecule 60 “single” and 30
“double” bonds are present, which can be modeled by extra electron charge density at the center of the double bonds compensated by electron charge deficiency at the single bond centers. Although similar, but considerably more complex, models have been proposed for CT0 [ 141 the present calculations indicate that including the van der Waals interactions alone is sufficient to correctly describe the orienta- . . . tional ordering m solid Cto. After the minimization the frequencies and eigenvectors of the lattice modes are calculated. The librational (f!) and translational (t ) inter-molecular modes in P2,/m CT,, are classified as I% 2% + 4B;+2B:t4A: and P=4%t+B;t2B:tA:. Of these 2 1 modes only the 12 gerade modes are Raman
active [ 4 1, whereas the IR spectrum will originate from the 9 ungerade modes. The calculated frequencies and the symmetry assignment for these modes are given in table 2 for the various potentials used

It is observed that the librational modes span the 15- 25 cm-’ range and that the translational modes are separated from these, spanning the 25-60 cm-’ range. In the first column of table 2 the experimentally determined frequencies are given. Although an assignment of these modes in the low frequency part (the “crowded” part) of the spectrum is not unambiguous, the overall agreement between the theoretical predictions and the experimental observations is good. Again the LJ ( 1) potential gives the best

fitting results. This is explicitly demonstrated in fig2 where the experimental Raman spectrum and the stick spectrum deduced from these data, as reported by van Loosdrecht et al. 141, are compared to the theoretical Raman spectrum obtained with the LJ ( 1)potential (see table 1). The predicted IR spectrum

calculated from the same LJ ( 1) potential is also shown in fig. 2, and awaits experimental conlirmation.

Conclusions
Lattice dynamics calculations on the low temperature structure of solid CT0 have been performed, using simple atom-atom van der Waals potentials to describe the interactions between the rigid CT0 molecules. The theoretically obtained crystal structure agrees with the experimentally determined low temperature structure; the symmetry is found to be P2,/ m. The calculations yield more insight in the orientational ordering, which involves a fixed rotation
of the four molecules in the unit cell around their long axes, and which can be described by a single angle with 96% 6.5 ‘. The calculated frequencies of the 12 Raman active lattice modes agree quantitatively with the recently observed Raman spectrum. The frequencies of the nine complementary IR active lattice modes are predicted

Lattice dynamics of KI perturbed by Jahn-Teller centres

Lattice dynamics of KI perturbed by Jahn-Teller centres

The theoretical study of the physical properties of JT centres (impurities, defects, etc. )embedded in a crystal is usually done in the framework of the cluster model [I]. In such a model the localized electron of the centre, which is in a degenerate state by definition of the JT effect, is assumed to interact only with a single localized vibrational mode, any lattice dispersion of frequency being neglected. This model is known to take into account several important features of the JT induced properties [I] but not the details of the optical spectra (Raman, infrared, vibronic sidebands of the absorption/emission spectra), which are very sensitive to the dynamics of the host crystal.

In a recent paper [2] the problem of the lattice dynamics perturbed by a JT centre was studied and resolved in the weak interaction limit by using the Green function method. It was found that a JT centre generates a very peculiar perturbation on the dynamical matrix.

On the other hand, a similar result can be obtained when the absorption lineshape due to a singletmultiplet transition (the dynamic JT effect) is studied. In this case the JT perturbed lattice dynamics is that relative to the excited electronic degenerate state. In both cases we conclude that the JT interaction induces a perturbation on the force constant matrix, which can be eventually seen in the structures of the phonon densities. Such a result, which will be presented in details in a forthcoming paper, can be deduced in the limit of the weak interaction from the absorption lineshape theory proposed some years ago (ref. [3] hereafter called MT).
In the present paper, after some consideration on the previous point, we present the computed one-phonon spectra for a (e x E) JT centre. As usual, (e x E) means that the electron is in a e-symmetry degenerate state and interacts with the E-symmetry displacements of the neighbour ions. Our aim is to see to what extent the JT dynamical perturbation can qualitatively modify the host lattice dynamics. So we do not apply it to any real impurity.
In the following we assume : 1) The JT centre is embedded in an ionic crystal. 2) The JT interaction is weak. 3) Only the n.n. of the JT centre are involved in the JT interaction.

The shape of the optical spectra is related to the symmetrized (E-symmetry in this case) projected phonon density p,(w) relative to the lattice dynamics, either of the ground electronic state (Raman and infrared spectra), or of the excited electronic state (vibronic sidebands). We remember that p,(o) is given by the imaginary part of the E-symmetry phonon propagator DE(a). So we first show the relevant JT induced processes which dress the unperturbed onephonon
propagator, and then we evaluated p,(o) from D,(o). By unperturbed propagator we mean the onephonon propagator when all the dynamical perturbations induced by the centre (change of mass and of force constants) are considered, but the JT interaction. It follows that the usual separation of the total Hamiltonian H in three terms (JT electron He, perturbed lattice HL and their interaction He,) can be used :

1) The electronic JT eigenfunctions and eigenvalues are evaluated in the static lattice and the electron is assumed to be in a doubly degenerate stationary state, whose energy E, falls in the host crystal energy gap. Such a state is either the ground state or the excited state, depending on the problem.

2) HL is the lattice Hamiltonian in the harmonic approximation, given by

Here I labels the ions in the crystal (I = 0 indicates the position of theJT centre) ; A,,. = A;. + 6A,,. is the force constant matrix of the perturbed lattice (the perturbation 6A,,. is zero for anisotopic centre); M, are the masses ; u, and p1 are the canonical dynamical variables (displacement u, and momentum p,). A labels the normal modes of the imperfect crystal whose frequency is on ; b, and b: are the annihilation and creation operators of a phonon on the A mode.

3) HeL is the JT interaction Hamiltonian. He, is usually written as follows

where u,, and u,, are the two E-symmetry combinations of the n.n. displacements of the JT centre; h(E1) = ox and h(E2) = - o,, where ox and o, are the Pauli matrices [I]. c is the coupling constant. HeL is very peculiar : it is localized both in the wavefunction space (He, # 0 only inside the degenerate e-space) and in the lattice space (only the n.n. of the JT centre interact with the electron). The complexity of the JT problem is immediately shown by evaluating the commutators of the electron and phonon operators entering He,, because :

i) He, does no commute with TL, because of the phonon operators ; [TL,H e =\ 0] This is the term which generates phonons (the so-called many-phonon process) also for normal impurities.

ii) The coupling coefficients h(E, i) are not commuting matrices : [h(E, 1), h(E,2) =\ 0]

iii) For the previous reason also

In the cluster model, where one assumes that only the n.n. of the centre are involved in HL, as they are in He,, only two degenerate E-symmetry modes with the same frequency are representative of the whole lattice dynamics. So the commutator (6) loses its meaning. From the properties of the commutators (4) and (5) one can see that, while the electron and the
phonons are truly independent when He, = 0, they are so entangled for HeL # 0 to determined new coupled states, the so-called vibronic states. When the vibronic states are studied in the cluster model, the role played by HeL on the electronic states is emphasized but the possible modifications of HL are neglected [I].

The opposite approach is here follow ed. In fact we compute the modifications of the one-phonon projected density of states pE(02), when the electron is supposed to be degenerate but structureless. Without the JT interaction, the unperturbed phonon propagator D;(t) (a dotted line in the diagram of figure 1) is given at T = 0 K by

where p;(02) is the unperturbed projected one-phonon density. (By &(w2) we mean its Hilbert transform.) In the harmonic approximation a phonon does not change frequency when a linear electron-phonon interaction is switched on and the only effect is a generation of many unperturbed phonons. When the electron is in a degenerate state, since the usual separation of the phonons from the electron (adiabatic approximation) is not valid anymore, the linear JT interaction modifies the frequency of the phonons. By following MT [3], one can prove that in the diagrams accounting for the many-phonon process there is always a class of higher order graphs (the crossing graphs) which give a contribution to the lower order graphs, when the electron line is removed by assuming that the electron is structureless. We draw such graph in figure 1 and decompose them as it follows :

In the previous graphs the squares are the commutators (6), the full line is the electron propagator and the dotted lines are the phonon propagators (7). The series on the last row can be summed in the weak interaction limit and the dressed E-symmetry propagator DE(m2) can be evaluated [2]. Note that, owing to the matrix character of W,, in the e-space, also D,(02) is a 2 x 2 matrix, so one finds that the resulting JT perturbed projected density pE(w) is composed
by two parts ; &o), each corresponding to one of the values of the frequency dependent JT dynamical perturbation 6A, (m)

to examine both points and emphasize the influence of 6A, on p$(o), we computed p,'(w) also for values of Q,, outside the range of validity of the very weak interaction limit.
The difference in the values here found for the average frequency on the distributions p:(o) corresponds to the splitting of E7 (o = 15 f 2 QJT), which is found in the cluster model in the very weak interaction limit. The computed pg(o) confirm that the more intense the JT perturbation, the greater the difference between ~ : ( oan)d p, (o).

As regard to the structure of p$(o), we note that &A, causes two effects. 1) It can generate a low frequency resonance whose intensity and frequency goes to zero for 611, 4 0. No gap modes are found in KI. 2) It strongly modifies the host lattice frequency distribution pg(o), mainly in the acoustic region. ,Therefore, we can suggest on a qualitative basis that 'the structures found in the one-phonon part of the experimental optical spectra of the ST centres, should correspond either to JT resonances or to peaks in the host lattice density of states perturbed by dA,. Moreover, since such peaks are not predicted in the cluster model, one must be careful in using it for interpreting the optical spectra.

The specific energy of the crystal deformed in accord with the matrix kε 1 was calculated as a function of the strain magnitude γ and was used in an equation that can be written as

Results and Discussion
3.1. Phase Transition. We obtained the static energy-volume (E-V) curves for the bcc, fcc, and hcp structures. From the static energy differences shown in Figure 1(a), we can see that the hcp structure is unstable in all ranges of volumes. The bcc structure is stable above 8 Å3. It means that the bcc Mo transforms to the fcc Mo under compression. By fitting the E-V data to the fourth-order finite strain EOS,39,40 we can obtain the enthalpy as a function of pressure. The enthalpy differences are shown in Figure 1(b). We find that Mo is stable in the bcc structure up to 703 ( 19 GPa, and then it transforms to the fcc phase. The hcp structure can not be the stable phase in the whole range of pressures. The equilibrium volume V0, bulk modulus B0, and its pressure derivatives B′ and B′′ are listed in Table 1. The agreement of our results with experimental data4,41 is very good. From Table 1, we can see that the bulk modulus of fcc
Mo at zero pressure is smaller than that of bcc Mo. It seems incompatible because the fcc phase is more densely packed than the bcc phase. Actually, the fcc phase of Mo is not observed at zero pressure and temperature. Comparing the first and second pressure derivatives of bulk modulus of bcc Mo with that of the fcc Mo, it is obvious that the bulk modulus of fcc Mo is larger than that of bcc Mo under ultrahigh pressure. Whether the bcc Mo is stable or not under high temperature and high pressure remains disputable. By calculating the Gibbs free energies of the bcc and fcc Mo from 350 to 850 GPa at room temperature up to 7500 K, Belonoshko et al. found that the bcc Mo transformed to the fcc phase below 700 GPa at high temperature.10 We also calculated the free energy within QHA and obtained the phase diagram of Mo at high pressure and temperature, which is shown as Figure 1(c). It can be seen that our bcc-fcc boundary above 700 GPa is consistent with that obtained by Belonoshko et al. Below 700 GPa, the bcc-fcc boundary from Belonoshko et al. only exists at the pressure range 350-700 GPa. However, from our calculations, the boundary can extend to 232 GPa (4000 K). It is very close to the SW data 210 GPa (4100 K), where the discontinuity in longitudinal sound velocity was detected.3 It may imply that the bcc Mo transforms to the fcc phase at this shock compression. If we extrapolate the bcc-fcc boundary to 119 GPa, the temperature will locate at 3557 K, a little larger than the latest double-sided laser-heated DAC melting temperature 3302 ( 140 K17 at this pressure. So we consider that the DAC melting curve might not be a simple solid-liquid phase boundary but the indication of an unusual transition. Indeed Santamarı´a-Pe´rez et al. also suggested that the melted Mo in DAC experiments was a complex structured liquid (may be tetrahedral and icosahedral structures with short-range order), more akin to the supercooled liquid phase than a simple liquid.17 Cazorla et al. also calculated the free energy of the bcc, fcc, and hcp structures and found that the hcp was the most stable phase at this high pressure and temperature region.11 However, their bcc-hcp boundary is very close to the present (and Belonoshko et al.) bcc-fcc boundary.

3.2. Elastic Properties. For each volume of the unit cell, the elastic constants (C11, C12, C44) of the bcc Mo were deduced from a polynomial fit of the strain energy for particular deformations listed in eqs 5-7. Ten symmetric values of γ in the range ( 5% were used to make the strain energy fit at each strain type. For P, the value was calculated from the static equation of state. The obtained elastic constants C11, C12, and C44 at equilibrium structure parameters are 472.3, 160.4, and 106.0 GPa, respectively. The results of these calculations under pressure are presented in Figure 2, which are in good agreement with the energy-dispersive X-ray diffraction data.42 But Duffy

et al. only measured the elastic constants up to 24 GPa, and they showed the pressure derivatives of elastic moduli (linear fit) were 7.3, 3.3, and 0.5.42 From Figure 2, we can see that the pressure dependence of elastic moduli under high pressure can not be described as linear relations. All the elastic moduli can be described as third polynomials perfectly, which are shown
as follows

C11 ) 472.3 + 6.06564P - 5.67 × 10-3P2 +
6.10951 × 10-6P3 (8)
C12 ) 160.4 + 2.71872P - 2.54 × 10-3P2 +
2.58221 × 10-6P3 (9)
C44 ) 106.0 + 2.06781P - 2.13 × 10-4P2 +
2.25069 × 10-7P3 (10)

The theoretical polycrystalline elastic modulus can be determined from the independent elastic constants. The average isotropic shear modulus G and bulk modulus B of polycrystalline can be calculated according to Voigt-Reuss-Hill approximations. 43 Then the isotropically averaged aggregate velocities can be obtained as follows

VP ) [(B + 4/3G)/F]1/2 (11)
VB ) (B/F)1/2 (12)
VS ) (G/F)1/2 (13)

where VP, VB, and VS are the compressional, bulk, and shear sound velocities, respectively. The bulk and shear modulus B and G at high pressure are shown in Figure 3(a), and the aggregate velocities are shown in Figure 3(b). As pressure increases, the B, G, VP, VB, and VS increase monotonously.

3.3. Phonon Dispersions. Figure 4(a) shows the obtained dispersion curves of bcc Mo at zero pressure along several highsymmetry directions in the Brillouin zone (BZ) for both transverse acoustical (TA) and longitudinal acoustical (LA) branches, together with the experimental dispersion curves.26,28 It can be seen that our results agree excellently with the
experimental data. Along Γ-H, the theory does not predict the

overbending at about q ) 0.8 in the TA branch. At the H point, the theoretic value is a littl larger (∼5%). Along Γ-P-H, the agreement in shape as well as absolute energies is very good. Along Γ N, the calculation of TA [110]〈001〉 captured the small softening near the zone boundary successfully. Our phonon branches at zero pressure spread up to 277 cm-1 at the N point. Meanwhile, the dispersion curves at 17.2 GPa are plotted in Figure 4(b). The phonon branches at this pressure spread up to 322 cm-1 THz at the N point. One finds the excellent agreement of this dispersion curve with the results from inelastic X-ray scattering (IXS).28 We repeated the phonon calculations for the other 13 different volumes. Some of the dispersion curves are plotted in Figure 5(a), which shows the well-known phenomenon in solids; i.e., the phonon frequencies increase as volume decreases, except the values along the Γ-P-H direction and Γ-N in the lower TA branch. As pressure increases, the softening of TA [110]〈001〉 along Γ-N becomes more and more obvious. Under ultra compression (V ) 7.69 Å3), the frequencies along P-H in TA branches soften to imaginary frequencies, indicating a structural instability. Actually, from our static structure calculations, the bcc phase transforms to the fcc phase at about 703 GPa (about V ) 8.07 Å3). To understand the stability of Mo, we also calculated the phonon dispersions of the fcc Mo in 14 volumes. Some phonon

dispersions are shown in Figure 5(b). It is obvious that the fcc Mo is unstable at zero pressure (V ) 16.14 Å3), and with the decrease of volume, all the frequencies of LA and TA increase. The calculation predicted the stability of the fcc Mo by promoting the frequencies of the phonons along Γ to X and Γ to L symmetry lines from imaginary to real. Though the phonon dispersions indicate the stability of the fcc Mo at 9.68 Å3, the enthalpy at this compression is also insufficient to overcome the potential barrier of the bcc-fcc phase transition. From the stabilities of the bcc and fcc Mo, one notes that the phase transition might not be induced independently by the dynamic instability. However, the phase transition occurs before the TA modes go to zero frequency, and the mode softening behaviors are related to the particular mechanism which is responsible for the phase transition.

Conclusions
We employed the density functional perturbation theory to investigate the phase transition, elastic properties, lattice dynamical properties, thermal equation of state, and thermodynamic properties of the bcc structure Mo. By comparing the static energy-volume and enthalpy-volume curves for the bcc, fcc, and hcp structures, we found that Mo is stable in the bcc structure up to 703 ( 19 GPa and then transforms to the fcc structure. The hcp structure can not be the stable phase in the whole range of pressures. With the QHA, the bcc-fcc boundary under high temperature and pressure was obtained. The present results qualitatively agree with those reported by Belonoshko et al. but need the confirmation by further experiments. Meanwhile, to obtain the more accurate phase diagram, the anharmonic phonon-phonon interactions and electron-phonon interactions must be considered completely. The calculated elastic constants at low pressure agree well with the experimental data. Under high pressure, the elastic constants can not be described as linear relations, but the third-order polynomial can
describe the functions of elastic constants versus pressure satisfactorily. From the elastic constants, the bulk and shear moduli and the sound velocities are obtained successfully. The calculated phonon dispersion curves accord excellently with experiments. Under pressure, we captured a large softening along H-P in the TA branches. When the volume is compressed to 7.69 Å3, the frequencies along H-P in the TA branches soften to imaginary frequencies, indicating a structural instability, while with pressure increasing, the phonon calculations on the fcc Mo predicted the stability by promoting the frequencies along Γ to X and Γ to L symmetry lines from imaginary to real. We also investigated the thermal properties including thermal EOS, thermal pressure, volume thermal expansion, and Hugoniot properties. From thermal expansion coefficients and heat capacity, we found that the QHA is valid only up to about the melting point at zero pressure, but under pressure, the validity of QHA can be extended to much higher temperature. On one hand, the QHA includes part of anharmonicity by allowing phonon frequencies to vary with volume. On the other hand, the pressure can suppress part of the anharmonicity by strengthening the bondings among atoms and lowering the vibration of atoms. The Debye temperatures calculated from the elastic constant agree well with that obtained from the Debye approximation.

Lattice Dynamics and Thermodynamics of Molybdenum from First-Principles Calculations

Lattice Dynamics and Thermodynamics of Molybdenum from First-Principles Calculations

of the ultrahigh pressure scale. The equation of state (EOS) of Mo at high pressure is being used as a calibration standard to the ruby fluorescence in diamond anvil compression (DAC) experiments.1 Because of its important position in the field of material science and condensed matter science, Mo has attracted tremendous experimental and theoretical interest in its wide range of properties recently. At ambient condition, Mo is in body-centered-cubic (bcc) structure and melts at 2890 K.2 But what is the most stable phase of Mo under ultrahigh pressure? Experimentally, the shock wave (SW) acoustic velocity measurements showed that there was a sharp break on the Hugoniot curve at about 210 GPa (at a calculated temperature of 4100 K), which indicated that a solid-solid phase transition occurred prior to melting at 390 GPa (at a calculated temperature of 10 000 K).3 To compare with the SW data, Vohra and Ruoff investigated the static compression of Mo by energy-dispersive X-ray diffraction and found that the bcc Mo was stable up to 272 GPa at 300 K.4 The phase transition in shock compression at 210 GPa was not observed. By further X-ray diffraction investigation, Ruoff et al. showed that Mo was stable in a bcc structure up to at least 560 GPa at room temperature.5 Theoretically, Moriarty suggested that Mo was stable in the bcc structure up to 420 GPa, where it transformed to a hexagonal close-packed (hcp) structure and then at 620 GPa to a face-centered closepacked (fcc) structure.6 Later, Boettger7 and Christensen et al.8 showed that the hcp phase of Mo was not stable, and the bcc phase transformed directly to the fcc phase at 700 GPa.

Belonoshko et al. confirmed the results of Boettger and showed the transition pressure was 720 GPa at zero pressure.9 By calculating the Gibbs free energies of the bcc and fcc Mo in the pressure range from 350 to 850 GPa at room temperatures up to 7500 K, Belonoshko et al. found that Mo had lower free energy in the fcc phase than in the bcc phase at elevated temperatures.10 Shortly after these new results were reported, Cazorla et al. found that the hcp Mo was noticeably more stable above 350 GPa at high temperature by calculating the Helmholtz and Gibbs free energies of the bcc, fcc, and hcp phases.11 The other intriguing problem is the melting properties. For the transition metals, such as Mo, Ta, and W, there are enormous discrepancies in melting curves between laser-heated DAC12-17 and SW3,18 methods. As for Mo (as well as Ta and W), several thousand degrees of discrepancies exist in extrapolating from DAC pressures of around 100 GPa12-15,17 to SW pressure of 390 GPa.18 As is known that the overestimation of the melting temperature exists in SW experiments, Errandonea13 corrected the SW data by considering 30% superheating. The revised melting temperatures are located at 7700 ( 1500 K (390 GPa), also much larger than the melting temperature (just above 4000 K) at this pressure extrapolated from DAC experiments. Results from the empirical and phenomenological melting models are dependent on the selection of the parameters.13,19-22 It is shown that the choice of different sets of parameters leads to huge differences in the melting temperatures at high pressure. In addition, theoretical results are consistent with the SW data at high pressure but diverge from DAC data below 100 GPa. All these results are inadequate to explain the extreme discrepancies in extrapolating the DAC data to the SW data.9-11,23 The highpressure melting curve of Mo still remains inconclusive up to now.

The low-temperature phonon spectrum of Mo as measured in inelastic neutron scattering experiments exhibited a variety of anomalies: the large softening near the H point, the T2 branch at the N point, and the longitudinal branch did

not display the rounded dip near q ) 2/3 [111], which was typical for “regular” monatomic bcc metals.24-26 With temperature increasing, the H point phonon displayed anomalous stiffening, which had been proposed to arise from either intrinsic anharmonicity of the interatomic potential or electron-phonon coupling.26 Theoretically, using the molecular dynamics (MD) simulations with environment-dependent tight-binding parametrization, Haas et al. reproduced the weakening of the phonon anomalies as the temperature increased.27 Later, Farber et al. determined the lattice dynamics of Mo at high pressure to 37 GPa using high-resolution inelastic X-ray scattering (IXS).28 Meanwhile, they calculated the quasiharmonic phonon spectrum up to the highest experimental pressure by linear response theory. Both the experimental and theoretical results showed an obvious decrease in the relative magnitude of the H point phonon anomaly under compression. Recently, Cazorla et al. obtained the phonon dispersion curves of the bcc Mo at zero pressure using the small displacement method,29 but under larger compression, there are no experimental and theoretical studies. The first-principles density functional theory is very successful in predicting the high-pressure behavior of phonon dispersion relations and their concomitant anharmonic effects.30 It is more effective to connect the phonon properties directly to the lattice dynamics under pressure, temperature, or their combination. One of the main purposes of this work is to investigate the lattice dynamics and thermodynamics of Mo under high pressure and temperature.

The aim of the present work is multiple. First, density functional theory with the generalized gradient approximation (GGA) has been used for first-principles studies of the phase transition and elastic properties of Mo under high pressure, and then the quasiharmonic approximation (QHA) has been applied to the study of the lattice dynamic properties, the thermal EOS, and the thermodynamic properties. The organization of this paper is as follows. In Section 2, we give a brief description of the theoretical computational methods. The results and detailed discussions are presented in Section 3. A short conclusion is drawn in the last section.

Computational Methodology

For many metals and alloys, the Helmholtz free energy F can be accurately separated as
F(V, T) = Estatic(V) + Fphon(V, T) + Felec(V, T) (1)

where Estatic(V) is the energy of a static lattice at zero temperature; Felec(V,T) is the thermal free energy arising from electronic excitations; and Fphon(V,T) is the phonon contribution. Both Estatic(V) and Felec(V,T) can be obtained from first-principles calculations directly. Density functional perturbation theory (DFPT) is a well-established method for calculating the vibrational properties from first-principles in the framework of QHA.31,32 Including part of the anharmonic effects by considering the volume dependence of phonon frequencies gives access to the thermal expansion, thermal EOS, and thermodynamic properties. Our calculations were performed within the GGA, as implemented in the QUANTUM-ESPRESSO package.33 A nonlinear core correction to the exchange-correlation energy function was introduced to generate Vanderbilt ultrasoft pseudopotential for Mo with the valence electrons configuration of 4d55s1. In addition, the pseudopotential was generated with a scalar-relativistic calculation using GGA according to the recipe of Perdew-Wang 91.34

During our calculations, we made careful tests on k and q grids, the kinetic energy cutoff, and many other parameters to guarantee phonon frequencies and free energies to be well converged. Dynamical matrices were computed at 29 wave (q) vectors using an 8 × 8 × 8 q grid in the irreducible wedge of the Brillouin zone. The kinetic energy cutoff Ecutoff was 60 Ry,
and the k grids used in total energy and phonon calculations were 20 × 20 × 20 and 14 × 14 × 14 Monkhorst-Pack (MP) meshes.35 The self-consistent calculation was terminated when the total energy difference in two successive loops was less than 10-12 Ry. A Fermi-Dirac smearing width of 0.02 Ry was applied for Brillouin zone integrations in phonon frequency calculations, and in the calculations of static energy and thermal electronic excitations we treated the smearing width of Fermi-Dirac as the physical temperature of electrons. The geometric mean phonon frequency is defined by

where ωqj is the phonon frequency of branch j at wave vector q, and Nqj is the number of branches times the total number of q points in the sum. With the tested parameters, the geometric mean phonon frequency was converged to 1 cm-1. The elastic constants are defined by means of a Taylor expansion of the total energy, E (V, δ), for the system with respect to a small strain δ of the lattice primitive cell volume V. The energy of a strained system is expressed as follows

where E(V0, 0) is the energy of the unstrained system with equilibrium volume V0; τi is an element in the stress tensor; and i is a factor to take care of the Voigt index.36 The independent elastic constants in cubic structure are C11, C12, and C44. For the calculations of three elastic constants, we considered three independent volume-nonconserving strains