Thermodynamic properties of wadsleyite with anharmonic effect

The thermodynamic properties of crystals can be routinely calculated by density functional theory calculations combining with quasi-harmonic approximation. Based on the method developed recently by Wu and Wentzcovitch (Phys Rev B 79:104304, 2009) and Wu (Phys Rev B 81:172301, 2010), we are able to further ab initio include anharmonic effect on thermodynamic properties of crystals by one additional canonical ensemble with numbers of particle, volume and temperature fixed (NVT) molecular dynamic simulations. Our study indicates that phonon–phonon interaction causes the renormalized phonon frequencies of wadsleyite decrease with temperature. This is consistent with the Raman experimental observation. The anharmonic free energy of wadsleyite is negative and its heat capacity at constant pressure can exceed the Dulong–Petit limit at high temperature. The anharmonicity still significantly affects thermodynamic properties of wadsleyite at pressure and temperature conditions correspond to the transition zone.


Introduction
Olivine is dominant mineral in upper mantle with volume fraction is about 61 %. It has two high-pressure phases: wadsleyite and ringwoodite. The transformation from olivine to wadsleyite, which is responsible to 410 km seismic discontinuity, and disassociation of ringwoodite to magnesiowüstite and perovskite, which causes 670 km seismic discontinuity, marks the begin and end of transition zone, respectively. 520 km seismic discontinuity is thought to result from the transformation from wadsleyite to ringwoodite. The metastable olivine phase transformations have also been regarded as one of main possible mechanisms of deep earthquakes (Kirby et al. 1996). The wedgeshaped zones of metastable peridotite probably persist to the bottom of transition zone or uppermost of lower mantle because of the transitions of olivine to its high pressure polymorphs kinetically hindered in cold slabs. Therefore, the basic properties of the olivine and its high-pressure phases and their phase relations are very critical to improve our understanding to upper mantle and deep earthquake and have been extensively studied (Akaogi et al. 1989;Bell and Rossman 1992;Fei et al. 2004;Gillet et al. 1991;Irifune et al. 1998;Katsura and Ito 1989;Li et al. 2007;Reynard et al. 1996;Shim et al. 2001;Weidner et al. 1984;Wu and Wentzcovitch 2007;Yu et al. 2007;Yu and Wentzcovitch 2006).
Density functional theory (DFT) calculations combining with quasi-harmonic approximation (QHA) have been widely used to calculate the thermodynamic properties of materials (Baroni et al. 2010;Wentzcovitch et al. 2010;Wu et al. 2008). QHA assumes that the phonon frequency only depends on the volume and temperature effect on phonon frequency comes from the volume variation. Because DFT calculations can precisely predict the phonon density of state of materials and QHA considers the volume effect on phonon frequency, QHA works excellent for MgO at temperature below 1,000 K at 0 GPa where intrinsic anharmonic effect can be neglected. However, intrinsic anharmonic effect, which results from phonon-phonon interaction and leads to temperature dependence of phonon frequencies at fixed volume, increases with the temperature. Above a certain temperature, the intrinsic anharmonic effect can dramatically modify the thermodynamic properties of material. For example, the thermal expansion of MgO and diamond from QHA calculations begin to deviate rapid the experiment result at about 1000 K (Karki et al. 2000;Pavone et al. 1993). This deviation is caused by intrinsic anharmonic effect ignored by QHA. Since the temperatures of the Earth's interior are usually far larger than 1,000 K, the effect of the intrinsic anharmonic effects on the thermodynamic properties and phase boundary of minerals need to be considered carefully.
In order to include the anharmonic free energy, usual methods need introduce several volume-sensitive parameters (Wu and Wentzcovitch 2009). These parameters are hard to determine. It is not easy to combine the method smoothly with DFT calculations. Recently, Wu and Wentzcovith (2009) developed a method that required only a constant to be determined. For materials with reliable high-temperature thermodynamic data at ambient pressure available, this constant can be easily determined by comparing thermodynamic properties from QHA and the experiment data. The dramatic improvement on the thermodynamic properties of MgO, diamond and forsterite at high temperature indicates that the method captures the anharmonic effect very well. The anharmonic free energy can be negative and positive depending on the materials.
The anharmonicity of forsterite has been discussed widely. Its anharmonic free energy is negative and causes that the heat capacity of forsterite at constant volume can cross the Dulong-Petit limit (Anderson et al. 1991;Bouhifd et al. 1996;Cynn et al. 1996;Gillet et al. 1991). In contrast, the anharmonicity of wadsleyite is not clear. The data of heat capacity at high temperature ([1,000 K), which are very useful to evaluate the anharmonicity, are not available for wadsleyite because of phase transition (Reynard et al. 1996). The only thermal expansion experimental data (Suzuki and Kumazawa 1980) indicates that QHA calculations overestimate the thermal expansion of wadsleyite. This seems to suggest that wadsleyite has positive anharmonic free energy, while wadsleyite and forsterite have different anharmonic features. But this conclusion is not robust. The experimental thermal expansion data are not very reliable in determining the anharmonic effect. For example, anharmonic free energy of forsterite can be positive or negative depending on which experimental thermal expansion data is adopted. Furthermore, Reynard et al. (1996) found that all measured phonon frequencies of wadsleyite decreased with temperature in fixed volume. This result indicates that the wadsleyite has negative anharmonic free energy, which is different from what got from thermal expansion. No reliable high-temperature thermodynamic data of wadsleyite are available to determine the unknown constant introduced in Wu and Wentzcovitch's method (Wu and Wentzcovitch 2009). This usually is true for high-pressure phases of other minerals. In order to apply the method developed by Wu and Wentzcovitch (2009) to these high-pressure phases, Wu further developed a method (Wu 2010) to calculate the unknown constant without requiring any experimental data. The result on MgO and diamond indicates that the method work excellently. Here, we will use the method to discuss the thermodynamic properties of wadsleyite.

Method
The Helmolhtz free energy within QHA is expressed as: where, q is wave vector and j is mode index. U 0 (V) is the internal static energy obtained by a first principles variable cell shape type structural optimization under hydrostatic pressure. The second and third terms are the zero-point (F zp H ) and thermal energies (F th H ), respectively. In the statically constrained QHA (Carrier et al. 2007), the crystal structure and phonon frequencies are unique functions of the volume. Temperature effects on these quantities are accounted for through volumetric effects only.
The temperature dependence of phonon frequency, the consequence of the anharmonicity, is expressed in implicit way in Wu and Wentzcovitch's method (Wu and Wentzcovitch 2009): where x is the phonon frequency from the first principles calculations and X is the temperature-dependence phonon frequency used in quasi-harmonic formula of the free energy. The temperature dependence of X is described by the temperature dependence of V 0 at a fixed volume V: where C is a constant, V and V 0 are the predicted quasiharmonic volumes at high and zero temperature under the same quasi-harmonic pressure. V À V 0 is the volume expansion caused by temperature variation at the pressure P. The two features of the anharmonicity, the anharmonicty (I) increases with temperature, and (II) decreases with pressure, are embodied in the term ðV À V 0 Þ=V 0 . This is why the method works well while there is only one constant C (Wu and Wentzcovitch 2009) introduced. The anharmonic free energy F anh is easily found to be: The anharmonic free energy can also be calculated through the molecular dynamics simulation by the thermodynamic integration: where U 0 ðVÞ is the minimal potential energy, UðT 0 ; VÞ is the potential energy at temperature T 0 , U qh ðT 0 ; VÞ is the potential energy in QHA at temperature T 0 . The brackets refer to ensemble average. Taking derivative with respect to temperature and using high temperature relation for the harmonic oscillator, at the high temperature. A can be calculated from one NVT ensemble first-principles molecular dynamic (FPMD) simulation and further used to determine C through Eq.
(2)-(6). The results in MgO and diamond (Wu 2010) indicates that C determined this way is volume and temperature insensitive, which support our assumption that anharmonicity can be well described by the last term in Eq.
(3) and C is a constant. Thermodynamic properties of MgO and diamond are significantly improved after including the anharmonic free energy. Computations were performed using Quantum Espresso, a package based on DFT plane wave, and pseudopotential (Giannozzi et al. 2009). The quasi-harmonic calculations of wadsleyite, which the current calculations are based on, have been published earlier (Wu and Wentzcovitch 2007). The local density approximation (LDA), with which the crystal structure are well predicted (Table 1) was used in the calculations. The pseudopotential for magnesium was generated by the method of von Barth and Car, while that for oxygen was generated by the method of Troullier and Martins (1991). The plane wave cut-off energy is 70 Ry for wadsleyite. Brillouin zone summations over electronic states were performed over 2 k 9 2 k 9 2 k mesh with (1/2,1/2,1/2) shift from origin. Dynamical matrices were computed on 2 q 9 2 q 9 2 q mesh using density functional perturbation theory (DFPT) (Baroni et al. 2001) and then interpolated in a regular 9 q 9 9 q 9 9 q mesh to obtain the vibrational density of state. FPMD were performed in the 112 atoms supercell with the same pseudopotential, exchange-correlation functional, and cut-off energy as QHA calculations. Time step for FPMD is 1 fs. The total simulation time is 4.5 ps. The last 3 ps simulation is used to get the ensemble average. The temperature was controlled by Andersen's method (Andersen 1980).

Results and discussion
The calculating A at 265.75 Å 3 /cell and 2,000 K is 0.030. This corresponds to a C -0.17. The anharmonic free energy of wadsleyite thus is negative (Fig. 1). Negative C means that the phonon frequency in fixed volume decrease with temperature. The temperature dependence of phonon frequency in fixed volume is expressed as (Wu and Wentzcovitch 2009) where i is mode index and b ¼ aðV; TÞV 2 K T V 0 V 0 K 0 . K T and K 0 is isothermal bulk modulus at T and 0 K, respectively, and V 0 ,V 0 , and V are described by Eq. 3. Since C is constant and the mode Grüneisen parameters, c i , are almost temperature independent, the temperature dependence of a i is mainly determined by b. As we know (Wu and Wentzcovitch 2009) that b increases nonlinearly with temperature below 500 K because of its dependence on thermal expansion, aðV; TÞ, and becomes almost temperature independent at higher temperature. Given that C is negative, phonon frequencies of wadsleyite will decrease with temperature nonlinearly below 500 K and almost linearly at higher temperature. Negative C is in consistence with Raman experimented by Reynard et al. (1996), which find that a i for all modes they measured are negative. Since b is approximately 2.8 9 10 -5 /K and 0:5\c i \1:4, the calculated a i is at the range of -0.5 to -1.3 9 10 -5 /K. This is close to the experimental a i at the range of 0 to -6 9 10 -5 /K given that the approximation and uncertainty involved in the experiment. The nonlinear behavior below 500 K in a i ignored by Reynard et al. may also be a source that leads to the difference between the calculated and measured a i . Equation (8) implies that the modes with larger c i vary quicker with temperature. This is consistent with Raman measurement in forsterite and wadsleyite (Gillet et al. 1991;Reynard et al. 1996). Negative C also implies that the heat capacity at constant volume will finally exceed the Dulong-Petit limit at high temperature.
The anharmonic effect on thermodynamic properties of wadsleyite is shown in Fig. 2. The anharmonic effect is not significant at low temperature. The thermal expansion, which is the most sensitive properties to anharmonic effect, almost is invariant with the anharmonic effect below 500 K. However, the anharmonic effect on thermodynamic properties increases rapidly with temperature above 1,000 K. As listed in Table 2, the anharmonic effect increases the thermal expansion 16 % and C P 4.3 % at 1,700 K at 0 GPa. Although the anharmonic effect on thermodynamic properties decreases dramatically with pressure, it is still significant at relevant mantle condition. As listed in Table 2, at 1,700 K and 13.5 GPa, where forsterite is transformed into wadsleyite, the anharmonic effect can increase the thermal expansion 9 % and thermal Grünisen parameter 6 %. The available experimental heat capacity data are also shown in Fig. 2. They are not enough to show whether or not the anharmonic method improves thermodynamic properties of wadsleyite since these data are below 1,000 K, The high-temperature thermodynamic data above 1,000 K, which are very useful to evaluate the anharmonic effect, are not available because wadsleyite is transformed into forsterite above 1000 K at 0 GPa (Reynard et al. 1996). Thermal expansion data from Suzuki  (1980) seem to support a positive C. However, thermal expansion data are not reliable based on our experience in forsterite, where we can get positive or negative C depending on which experimental data are adopted.

Conclusions
In conclusion, we ab inito investigated the anharmonic effect on thermodynamic properties of wadsleyite. Wadsleyite has negative anharmonic free energy. The calculated renormalized phonon frequencies decrease with temperature, which is consistent with the Raman observation. The anharmonic effect increases the thermal expansion of wadsleyite at high temperature. The experimental values measured by Suzuki and Kumazawa (1980) are significantly smaller than the calculated thermal expansion. Noted that the precise measurement of thermal expansion is still challenge even for forsterite, which is stable phase at 0 GPa, our results suggest that the experiment by Suzuki and Kumazawa (1980) may underestimate the thermal expansion of wadsleyite. Although the anharmonicity decreases with pressure, its effects on thermodynamic properties are still significant at the pressure and temperature conditions of transition zone. For example, the anharmonic effect can increase the thermal expansion 9 % and thermal Grünisen parameter 6 % at 1,700 K and 13.5 GPa.