First-principles calculations of high-pressure physical properties anisotropy for magnesite

The first-principles calculations based on density functional theory with projector-augmented wave are used to study the anisotropy of elastic modulus, mechanical hardness, minimum thermal conductivity, acoustic velocity and thermal expansion of magnesite (MgCO3) under deep mantle pressure. The calculation results of the phase transition pressure, equation of state, elastic constants, elastic moduli, elastic wave velocities and thermal expansion coefficient are consistent with those determined experimentally. The research results show that the elastic moduli have strong anisotropy, the mechanical hardness gradually softens with increasing pressure, the conduction velocity of heat in the [100] direction is faster than that in the [001] direction, the plane wave velocity anisotropy first increases and then gradually decreases with increasing pressure, and the shear wave velocity anisotropy increases with the increase of pressure, the thermal expansion in the [100] direction is greater than that in the [001] direction. The research results are of great significance to people’s understanding of the high-pressure physical properties of carbonates in the deep mantle.

Single-crystal elastic constants. The elastic properties of the earth's minerals are crucial to understanding their internal properties, especially in terms of their chemical composition and the propagation of seismic acoustic waves. Magnesite has six ( c 11 , c 12 , c 13 , c 14 , c 33 , c 44 ) independent elastic constants since c 66 = (c 11 − c 12 ) 2 . In order to confirm its mechanical stability, the following mechanical stability criteria are checked 26 : In this work, all the calculated elastic stiffness constants c ij satisfy the mechanical stability criteria, so it may be said that magnesite is mechanically stable.
The calculated elastic constants of magnesite from 0 to 80 GPa are plotted in Fig. 3 and the data at 0 GPa are summarized in Table 2, compared with the previous experimental 15 and theoretical 11,16 results. It can be clearly seen from Fig. 3 and Table 2 that the present calculated elastic constants of magnesite are in excellent agreement with the available experimental and theoretical results, and gradually increase with the pressure.
Anisotropy of elastic modulus. The polycrystalline elastic moduli, such as bulk modulus B , shear modulus G and Young's modulus E , can be evaluated by Voigt-Reuss-Hill scheme [27][28][29] . For rhombohedral magnesite, (1) c 11 > |c 12 |, c 44 > 0, 2c 2 13 < c 33 (c 11 + c 12 ), 2c 2 14 < c 44 (c 11 − c 12 ).     According to the bulk modulus B and shear modulus G , Young's modulus E is defined as E = 9BG (3B + G). Figure 4 presents the changes of bulk modulus B , shear modulus G , and Young's modulus E of magnesite along with the previous experimental 15 and theoretical 11,16 results with pressure. As shown in figures, the present calculated elastic moduli increase smoothly and monotonically with increasing pressure, which agree well with the experimental and theoretical data.
The elastic anisotropy in mineral is of great significance due to its implication in geoscience as well as in crystal physics. In order to evaluate the elastic anisotropy of magnesite, Ranganathan and Ostoja-Starzewski universal anisotropy index 30 , Kube's log-Euclidean anisotropy index 31 , and Chung and Buessem percent elastic anisotropy 32 are used. The A U ,A L ,A B , and A G are given by the following relations: (2) B V = 2c 11 + c 33 + 2c 12 + 4c 13 9 (3) G V = (2c 11 + c 33 ) − (c 12 + 2c 13 ) + 3 2c 44 + (c 11 − c 12 ) 2 5 (4) B R = 1 (2s 11 + s 33 ) + 2(s 12 + 2s 13 ) (5) G R = 15 4(2s 11 + s 33 ) − 4(s 12 + 2s 13 ) + 3(2s 44 + s 66 )  Fig. 5a may be observed that A U and A L increase with the increase of pressure, and the change trend is basically the same. It is found in Fig. 5b that the percentage of shear anisotropy increases with the increase of pressure, and the percentage of bulk anisotropy decreases, and the increase in the percentage of shear anisotropy is much greater than the decrease in the percentage of bulk anisotropy, this means that the contribution of shear anisotropy in the elastic anisotropy of magnesite is greater than that of bulk anisotropy.
In order to furthermore elucidate this anisotropic behavior, the most straightforward method is to plot the three-dimensional contours of mechanical moduli. The direction dependent shear modulus ( G ) and Young's modulus ( E ) for rhombohedral crystals can be defined as: where s ij are the usual elastic compliance constants and l 1 , l 2 , and l 3 are the direction cosines in any arbitrary direction. The ElasticPOST program 33,34 is used to obtain the 3D spatial distribution and their projection of shear modulus and Young's modulus for magnesite at various pressures, and the results are displayed in Figs. 6 and 7, respectively. As can be seen, the 3D figures of shear modulus and Young's modulus reveals a large degree of deviation in shape from the sphere. This means that magnesite has a strong anisotropy, which also confirms the calculation results in Fig. 5. The comparative analysis of shear modulus and Young's modulus for different directions as seen from the planar projections also indicates the anisotropy level.
Anisotropy of mechanical hardness. Vickers hardness is a fundamental property that is essential to describe the mechanical behavior of mineral, various semi-empirical relations have been proposed to estimate hardness using the elastic moduli. Vickers hardness is predicted using two theoretical models of hardness: Chen's model 35 : The calculated Vickers hardness of magnesite are depicted in Fig. 8 from 0 to 80 GPa. As illustrated in Fig. 8, the Vickers hardness decrease with increasing pressure, indicating magnesite becomes softer under high pressure, the Vickers hardness predicted by the Chen's model is smaller than that Tian's of the model in the entire pressure range. In order to evaluate the anisotropy of Vickers hardness of magnesite. The direction dependent hardness ( H ) can be obtained by fitting the direction dependent bulk modulus ( B ) and Young's modulus ( E ), defined as: H = 0.130548175274347E 2.2484942942017 B −1.51675853808829 , where B = 1 (s 11 + s 12 + s 13 ) − (s 11 + s 12 − s 13 − s 33 )l 2 3 . The 3D spatial distribution and its projection of Vickers hardness for magnesite at various pressures are presented in Fig. 9. The Vickers hardness exhibit strong direction-dependent changes, resulting in large anisotropy, The 2D representations planar projection in different directions also show this result. The elastic wave velocities of magnesite are shown in Fig. 10 from 0 to 80 GPa. Figure 10 show that the calculated elastic wave velocity is in good agreement with the previous experimental results 15 within the studied pressure range, the plane wave velocity v P propagate more speedily than the shear wave velocity v S . The consistency between the calculated elastic wave velocity and the experimental results provides reliability for further research on elastic wave velocities anisotropy.
Directional elastic wave velocities are computed by solving Christoffel's equation det C ijkl n j n l − ρv 2 δ ik = 0 38 , where C ijkl are the elastic stiffnesses, the n j are unit vectors of the wave propagation direction, v is the acoustic velocity, and δ ik is the Kronecker δ . Using AWESoMe program 39,40 with quadruple precision, the plane wave velocities and shear wave velocity and the shear wave splitting of magnesite in different propagation directions under various pressures are obtained, 3D representation of the elastic wave velocity and the shear wave splitting of magnesite are plotted in Fig. 11. It is observed from Fig. 11(left) that the plane wave velocities have minimum values along the z direction, firstly decreasing with the increase of pressure, and then gradually increasing. For www.nature.com/scientificreports/ the two shear wave velocity (fast and slow), the minimum values of the two wave velocities are shifted from the z direction yet they are still allocated around the z direction, but the magnitude of the shift gradually increases with the increase of pressure, especially the fast shear wave velocity, the results of shear wave splitting in Fig. 11(right) also further verify this result. Anisotropy of plane and shear wave velocity can be defined as A P = v P,max − v P,min v P × 100% and A S = v S,max − v S,min v S × 100% , respectively. The calculated elastic wave velocities anisotropy of magnesite is presented in Fig. 12 and the data at 0 GPa are listed in Table 3, along with the previous experimental 15,41,42 and theoretical 1,11 results. It can be found from Table 3 that the maximum error between the calculated plane wave velocity anisotropy and the experimental 15 value at 0 GPa is about 2.5%, and the maximum error between the shear wave velocity anisotropy and the experimental value 42 is about 2.75%, indicating that the calculated data are in agreement with available experimental data. At low pressure, the plane wave velocity anisotropy increases with the increase in pressure, but gradually decreases at high pressure. However, the experimental result of Yang et al 15 is that the plane wave velocity anisotropy increases with increasing pressure. The shear wave velocity anisotropy Anisotropy of minimum thermal conductivity. The thermal conductivity is a measure of material's heat conduction ability. Generally, the thermal conductivity decreases to a limit value considered as the minimum thermal conductivity with increasing temperature. Therefore, it is of great significance to study the minimum thermal conductivity of magnesite. The minimum thermal conductivity of magnesite is calculated on the basis of Clark's model 43 and Cahill's model 44 . In the Clarke model, the minimum thermal conductivity can be thought of as the limit the average phonon mean free path → the interatomic spacing. The Cahill model instead use a wavelength dependentmean free path to incorporate wave mechanics in the description of the average phonon mean free path. These models work well for many materials and give an intuitive description of the phonon limit of thermal conductivity. Clark's Model: Cahill's Model: where M a is the average mass per atom, E is Young's modulus, ρ is the density, M is the molar mass, n is the atomic number density per unit volume, k B is Boltzmann's constant, N A is Avogadro's number, respectively. Based on the two theoretical model, the calculated minimum thermal conductivity of magnesite from 0 to 80 GPa is shown in Fig. 13. It is seen that the minimum thermal conductivity of magnesite increases with the increase of the external  To investigate the anisotropy of thermal conductivity, which can be summed from the plane wave velocities ( v P ) and two shear wave velocity ( v S1 and v S2 ). Therefore, the expression of Cahill's model can be changed in form as follows: The calculated minimum thermal conductivities of magnesite in principal directions are also presented in Anisotropy of thermal expansion. The thermal expansion coefficients and their temperature-pressure dependence are of importance in estimating the thermal properties of minerals. In present work, The Debye quasi-harmonic approximation (QHA) is used to calculate the thermal expansion coefficients of magnesite 45 . The volumetric thermal expansion coefficient ( α V ) can be obtained by the following expressions: where γ , C V , B T , V and θ D represent the thermal Grüneisen parameter, the heat capacities, the isothermal bulk modulus, the volume and the Debye temperature, respectively. The volume thermal expansion coefficient of magnesite at 300 K and 0 GPa is 3.376 × 10 -5 /K, in good agreement with the present calculated value of 3.688 × 10 -5 /K. Having obtained the volumetric thermal expansion at different temperatures and pressures, the thermal expansion along different directions can be calculated from the linear compressibility. For rhombohedral crystal, the expressions are as follows 46  The anisotropic linear thermal expansion coefficients of magnesite at various pressures are calculated and are depicted in Fig. 14. As can be seen, the thermal expansion in the [100] direction is the largest relative to the [001] directions in magnesite, and it decrease with increasing pressure. Unfortunately, there is no experimental data or theoretical calculation results to compare with the linear thermal expansion coefficient of magnesite. Thus, the present work is beneficial for future research on the thermal properties of minerals.

Conclusions
The anisotropy of elastic modulus, mechanical hardness, minimum thermal conductivity, acoustic velocity and thermal expansion of magnesite under high pressure are investigated using the first-principles calculations within the density functional theory. The calculated phase transition pressure, equation of state, elastic constants, elastic moduli, elastic wave velocities and thermal expansion coefficient of magnesite are in excellent agreement with the previous experimental and theoretical results. It provides reliability for further research on the anisotropy of elastic modulus, mechanical hardness, minimum thermal conductivity, acoustic velocity and thermal expansion. The results of shear modulus and Young's modulus show that magnesite has strong anisotropy. The Vickers hardness changes strongly in different directions, leading to large anisotropy and softening under high pressure. Due to the higher probability of phonon collision in the [100] direction, the minimum thermal conductivity in the [100] direction is higher than that in the [001] direction and increases with the increase of pressure. The propagation of the plane wave along the z direction has a minimum value, which decreases first and then gradually increases as the pressure increases. The minimum value of the two shear wave velocities shifts from the z direction, and the magnitude of the shift gradually increases with the increase of pressure, especially in the fast S-mode. The plane wave velocity anisotropy first increases and then gradually decreases with increasing pressure, and the shear wave velocity anisotropy increases with the increase of pressure. As discussed in literature 1 , the elastic anisotropy of magnesite is much greater than that of the main minerals in the mantle, and its local enrichment provides a new explanation for the large local anisotropy in the transition zone. Finally, the anisotropy of thermal expansion is studied using the Debye quasi-harmonic approximation and elastic constants. It is found that the anisotropic linear thermal expansion coefficients in the [100] direction is the largest relative to the [001] directions and decrease with increasing pressure. The present work helps people to further understand the high-pressure physical properties of magnesite under deep mantle conditions, and also has important geophysical significance.

Methods
First-principles calculations based on density functional theory 48,49 are performed by using a Vienna Ab Initio Simulation Package (VASP) 50,51 with the projector-augmented wave method (PAW) 52 . The Perdew-Burke-Ernzerhof revised for solids(PBEsol) in GGA 53 is used to expound the exchange-correction function and calculate the self-consistent electronic density. The valence electron configurations are chosen 2p 6 3s 2 for Mg, 2s 2 2p 2 for C, and 2s 2 2p 4 for O. Based on the results of plane-wave cutoff energy and k-mesh convergence tests, the cutoff energy for plane wave extension of the R 3 c and C2/m structures for MgCO 3 are set to 850 eV and 880 eV, and the Brillouin zone of Monkhorst-Pack grid sampling 54 is 9 × 9 × 2 and 4 × 5 × 5, respectively. The convergence threshold for electronic self-consistent field and forces acting on the atoms are 1.0 × 10 −8 eV and 0.02 eV/Å, respectively. The elastic constant is obtained using the stress-strain method 55,56 . The thermodynamic properties are calculated by the quasi-harmonic approximation (QHA) Debye approach 45 .