Theoretical study of anisotropic structural, electronic, mechanical and thermodynamic properties of rare-earth (R = Y, La) oxysulfides

Computational Materials Science 125 (2016) 154–167

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Theoretical study of anisotropic structural, electronic, mechanical and thermodynamic properties of rare-earth (R = Y, La) oxysulfides Y.F. Li a,⇑, B. Xiao b, Y.M. Gao a, Y.H. Cheng c a

State Key Laboratory for Mechanical Behavior of Materials, Xi’an Jiaotong University, Xi’an 710049, China Department of Earth Sciences, University College London, London WC1E 6BT, England, United Kingdom c State Key Laboratory of Electric Insulation and Power Equipment, Xi’an Jiaotong University, Xi’an 710049, China b

a r t i c l e

i n f o

Article history: Received 4 June 2016 Received in revised form 1 August 2016 Accepted 29 August 2016 Available online 10 September 2016 Keywords: Rare-earth oxysulfides Anisotropy Electronic structure Thermal expansion

a b s t r a c t The anisotropic structural, electronic, mechanical and thermodynamic properties of rare-earth (R = Y, La) oxysulfides are calculated by first-principles using density functional theory. From the calculated electronic structure and electron density topology, we find that the crystal structure of R2O2S is dominated by strong polar RAO and RAS bonds, and the ionic bond is the main bonding mechanism. The negative formation enthalpies imply that the formation of R2O2S from reactants (R2O3 + H2S) is energetically favorable. The two oxysulfides show obvious anisotropy in elastic properties along different crystallographic directions, where Young’s modulus shows stronger anisotropy than bulk modulus. The anisotropy in sound velocity is obtained by solving the Christoffel equation, the sound modes are isotropic at the basal plane (x-y plane), and weak anisotropy is seen on (1 0 0) crystallographic plane. The temperature dependence of volumetric TECs of both oxysulfides are similar, and the linear TECs in [0 0 1] direction are slightly higher than that of the [1 0 0] direction. At room temperature, the computed average linear TECs for Y2O2S and La2O2S are 9.5
106 K1 and 9.4
106 K1, respectively. Using Slack’s model, it is found that thermal conductivities are 10.1 and 5.3 W/m K in a range from 300 K to 1500 K for Y2O2S and La2O2S, respectively. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction The rare-earth oxysulfides R2O2S (R = Y and La) have raised great research interests for phosphors, additive, catalyst and so on during the past decades. They are wide-gap semiconductors and have been used in a variety of applications, e.g., Y2O2S: Tb has used in X-ray imaging detectors [1]; and properly doped La2O2S is an excellent phosphorescent material for lasers [2]. Another important application is in steelmaking industry, the oxygen and sulfur elements always deteriorate the properties of steel, i.e., the sulfur element leads to the brittle of steel in hot-rolling temperature by forming low melting point FeS/Fe eutectic. Therefore, it is a big scientific issue for removing impurities (oxide and sulfide) in steelmaking practice. The high chemical activity of rare-earth elements and their strong affinity with oxygen and sulfur make them the most promising candidates of additives [3]. Actually, considerable quantity of R2O2S is formed by adding rare-earth elements in the liquid steel to inhibit grain growth as well as to reduce the oxide and sulfide impurities. As a result, the high tem⇑ Corresponding author. E-mail address: [email protected] (Y.F. Li). http://dx.doi.org/10.1016/j.commatsci.2016.08.050 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

perature mechanical and corrosion properties of the steel can be improved significantly [4]. Until recently, elastic, electronic and optical properties of R2O2S have been investigated. Li and Ahuja [5] reported Y2O2S is an indirect band gap material with the band gap value of 3.0 eV, and they also provided the bulk modulus of 125 GPa and dielectric constant of 5.3. Mikami and Oshiyama [6–8] claimed that the calculated band gap of Y2O2S (2.61 eV) is significantly smaller than experiment (4.8 eV); later, they also provide the calculated band gap of La2O2S (3.4 eV). Vali [9] found that La2O2S is also an indirect band gap semiconductor, and the origin of anomalously Born effective charges and the bonding character were also discussed. Itoh and Inabe [10] provided the reflectivity spectra of single crystals of Y2O2S at 8 and 300 K using synchrotron radiation as light source; a sharp exciton absorption peak observed at 6.33 eV was associated with the transition from the S 3p state to the Y 4d state. Polfus et al. [11] studied the enthalpies of dissociative absorption of H2S or H2O into anion vacancies in La2O2S, they found oxygen vacancies were energetically favorable over sulfide ions vacancies over a quite broad range of environmental conditions. Although the electronic, bonding, and optical properties of R2O2S have been extensively investigated from practical viewpoint

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

as function materials. To our knowledge, the information about the anisotropic mechanical property, Grüneisen parameters, temperature dependence thermal expansion and thermal conductivity of them are scarce. Those properties are important to further understand the intrinsic behaviors of R2O2S in cast irons and steels. In this work, we will perform systematic studies on the electronic structures, charge density topology, and elasticity of R2O2S crystal structures. Phonon spectra and other thermal physical properties derived from it are also calculated and discussed. 2. Methods and details 2.1. Density functional theory calculations The crystal structure of R2O2S is illustrated in Fig. 1a, which adopts trigonal structure (P3m1) and consists of a [R2O2] parallelogram basic plane with an angle of 75°. The basic plane is parallel to the [0 0 1] crystal direction and its normal vector is the [1 1 0] direction (Fig. 1b). From Fig. 1c, a periodic stacking of extended two-dimensional [RAO] basic layer parallel to the (0 0 1) crystal plane can be observed, which are separated by a hexagonal S atomic layer in R2O2S cell. The [R2O2] parallelogram basic plane sites on the side surface of 2D [RAO] basic layer, while the frontview surface including six distorted edges does not form a plane with the edges angle 113°. Take the basic unit of 2D [RAO] layer out and shown in Fig. 1d, the distorted hexagonal prism can be seen clearly, also we can see that the 2D [RAO] basic layer are formed by corner-sharing RO4 tetrahedral, and the R atoms are coordinated to one axial and three equatorial O atoms [12]. All first-principles calculations were performed by CASTEP code with plane wave basis and periodic boundary condition (PBC) [13– 15]. The Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm was applied to relax the structure where both lattice parameters and atomic fractional coordinates were optimized simultaneously. The pseudo-ionic core and valence electrons interactions were described by ultra-soft pseudopotentials (USPPs). The valence electrons of USPPs were: Y: 4d15s2, La: 5d16s2, O; 2s22p4, S: 3s23p4, respectively. The exchange-correlation density functional was treated within the local density approximation (LDA) using CAPZ approximation to compute the lattice structure, mechanical and phonon properties. As we know, LDA may underestimate the lattice constants. However, for many layered materials where the inter-layer interactions may play the important role for the equilibrium geometry, LDA could accidentally give reasonable equilibrium structural parameters [16]. Besides, it is well known that LDA may not reliable to describe the band structures (band gaps) in rare-earth compounds [17,18].

155

In this work, we employed the hybrid density functional sX-LDA for electronic structure calculations under the framework of generalized Kohn-Sham density functional theory (GKS-DFT). The exchange energy functional of sX-LDA [19–21] divides into shortand long-range parts and replaces LDA short-range exchange part totally by nonlocal Hartree-Fock exchange, screened using the Thomas-Fermi function. Meanwhile, the correlation energy functional is the same as that of LDA. In contrast to the bare ThomasFermi approximation, this method takes advantage of LDA for treating the electron-electron correlation interactions [19]. The kinetic energy cut-off value of 600 eV was used for plane wave expansion in reciprocal space. The sampling k-point mesh in the first irreducible Brillouin zone was generated by Monkhorst-Pack method [22], using a 8
8
4 grid for each structure. The convergence criterion for electronic self-consistency was 108 eV per atom for total energy, and the forces per atom were reduced to 1
104 eV/Å. The finite element displacement method was applied to calculate the phonon spectra by setting a cutoff radius of 3 Å. The dynamical matrices at given q points in first BZ can be obtained from the Fourier interpolation algorithm, the quality of the q factor spacing was set to 0.015 Å1, and the BZ path was set as C-A-H-K-C-M-L-H in this study. The phonon frequencies are solutions of the secular equation of the dynamic matrices at the q point [23].

2.2. Thermodynamic properties In order to calculate the thermal expansion coefficient (TEC), obtaining the Helmholtz free energy at the elevated temperature is essential, which can be calculated by F(V, T) = Egs(V) + Fvib(V, T) + Fele(V, T) + EZPE in general. The ground state total energy (Egs) is given by standard DFT total energy calculation at 0 K. The vibrational free energy (Fvib) is calculated by phonon spectra in quasiharmonic approximation (QHA). From the calculated phonon density of states (PHDOS) at the given volume V, the lattice vibration free energy Fvib can be obtained [24]. Fele is the thermal electronic contribution to Helmholtz free energy, which is obtained by Mermin statistics Fele = Eele  TSele. The thermal electronic contribution is neglected for rare-earth oxysulfides, because they are wide band gap semiconductors. Besides, the zero point energy (ZPE) can be obtained by integration of the phonon DOS [25]:

EZPE ¼

1 2

Z
hxgðV; xÞdx

ð1Þ

where gðV; xÞ is the PHDOS at each volume. Minimizing the Helmholtz free energy F (V, T) with respect to cell volume V at different

Fig. 1. Schematic diagrams of R2O2S (R = Y, La) crystal structures: (a) R2O2S crystal structure, (b) RAO basic plane, (c) Extended RAO layer, (d) Basic unit of extended RAO layer. The blue, pink, and yellow balls represent R, O, S atoms, respectively. The bond angles and bond lengths (in Å) outside the parenthesis are indicated for Y2O2S, while the values inside the parenthesis correspond to La2O2S. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

156

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

temperatures gives the equilibrium volumes. This can be fitted by Birch-Murnaghan isothermal equation of state:

" #  ð1B00 Þ Bo V o 1 V V B00 þ Eo EðVÞ ¼ 0 þ þ V o 1  B00 B0 B00  1 V o

ð2Þ

Finally, the equilibrium volume V(T) at T is obtained. The volu T metric TEC k(T) can be determined by k ¼ V10 @V , where V0 repre@T sents the equilibrium volume at 0 K. In order to evaluate the chemical stability of bulk R2O2S phases, we have calculated the formation enthalpy of them from the reaction R2O3 (S) + H2S (S) ? R2O2S (S) + H2O (S) by the following expression.

Dr H ¼ Etot ðR2 O2 SÞ þ Etot ðH2 OÞ  Etot ðR2 O3 Þ  Etot ðH2 SÞ

ð3Þ

Here, Dr H is the formation enthalpy. Etot(X) is the total energy of the X phase according to the standard DFT total energy calculation at 0 K; and X denotes the R2O2S, H2O, R2O3, or H2S phases. In this work, H2O, R2O3, and H2S adopt space groups P63/mmc, P3m1, and Pbcm, respectively. In our current work, most thermal physical properties and the analyzing of charge density topology were obtained using our home-made post processing programs. 3. Results and discussions 3.1. Stabilities and electronic structures The relaxed atomic positions and lattice parameters of both rare-earth oxysulfides are listed in Table 1. Comparably speaking, the average deviation of lattice parameters between experimental data and LDA results are less than 2%. LDA method provides accurate cell volumes, which are close to the experimental data. As we know, the accurate equilibrium volume is important for calculating thermal expansion. If the predicted cell volume is bigger than experimental value, the force constants which originate from chemical bonding may be underestimated, and thereby the phonon frequency probably would be smaller than experimental data. The formation enthalpy was calculated and shown in Table 1. Negative values of Dr H suggest that the formation of R2O2S from reactants (R2O3 + H2S) is energetically favorable. Polfus et al. [11] also claimed that the corresponding reaction is quite exothermic for the formation of La2O2S. Haynes and Brown [28] pointed out that La2O2S can be prepared by heat treating La2O2SO4 and a few percent H2 and H2S at 1273 K. Although the electronic bonding characters of rare-earth oxysulfides have been analyzed via first-principles calculations [5–9],

we still want to present here a brief discussions by means of structural motifs (Fig. 2), electronic structures (Fig. 3) and electron density maps (Fig. 4). Although, the fundamental structural blocks are illustrated in Fig. 1, and a short discussion has been given in previous section, it is still useful to analysis the local coordination number of each element in the structure from the point view of chemical bonding. In Fig. 2, the coordination species of R (Y, La), O and S is displayed. In R2O2S structure, R3+ cations have a coordination number seven with four O atoms and other three S atoms forming a capped trigonal antiprism. All three RAS bonds are symmetrically equivalent. Among four RAO bonds, three equatorial RAO bonds are identical, and the remaining apical RAO bond aligned in [0 0 1] crystallographic direction is longer than those equatorial bonds. Moreover, this apical RAO bond is also perpendicular to RS3O3 trigonal anti-prism. On the other hand, the coordination species of O atom is a distorted OR4 tetrahedral, and that for S atom is a distorted SR6 octahedral, as shown in Fig. 2b and c. In OR4 tetrahedral, the local symmetry of O atom is reduced from Td to C3v, because the RAO bond pointed in [0 0 1] direction is not the same as other three. In the case of S atom, all six S-R bonds are identical. However, RASAR angles do not match with those of ideal octahedral, but they are consistent with a trigonal antiprism with point group C3v. From the point of view of orbital hybridization, it is easy to argue that O and S atoms could form sp3 and sp3d2 type hybridized orbitals. Since Y and La atoms have complicated valence shells ((n  1)d, ns and np), i.e., the energies of (n  1)d, ns and np atomic orbitals are fairly close to each other. Therefore, Y or La could also form sufficient number of hybridized orbitals (spd5 or sd5 + pz) to support a large coordination number like seven in a capped trigonal anti-prism. The bond lengths of OAO, SAO and RAR pairs are larger than 3.0 Å, the strong interaction is not expected in the current structures. In order to discuss the chemical bonding mechanism in a qualitative way, we show the calculated band structures and densities of states of Y2O2S and La2O2S in Fig. 3, using sX-LDA exchangecorrelation functional. From the band structures of both oxysulfides, the top of the valence band is situated at the C points and the bottom of the conduction band at the K point. The band gap of 4.2 eV for both Y2O2S and La2O2S indicate they are semiconductors, which are close to the experimental data (4.8 eV) [6,27]. From the calculated angular-momentum projected density of states, the orbital hybridizations of O and S atoms are almost absent, i.e., their s and p states do not mixing and only form two separated bands. On the other hand, there is sufficient mixing of (n  1)d and ns states of R (Y or La) atom with O and S atoms. At the same time, we also find the strong mixing of (n  1)d and ns orbitals of Y or La atom. Therefore, the calculated electronic density of states implies the orbital hybridization occurs only in Y or La

Table 1 The calculated atomic fractional coordinates, lattice parameters (a, c, in Å; V0 in Å3; q in g/cm3), formation enthalpy (eV/f.u.) and band gap (eV) of R2O2S (R = Y, La). Compounds

Atomic positions1 u

Y2O2S La2O2S 1 a b c d e f g h i j

Lattice parameters v

a

0.282 (0.282 ) 0.278 (0.279a)

a a

0.630 (0.631 ) 0.629 (0.629a)

q

V0

a

b

3.73 (3.75 , 3.791 ) 4.04 (4.035a, 4.049f)

Energy parameters

Dr H

c a

b

6.48 (6.525 , 6.596 ) 6.90 (6.914a, 6.939f)

b

78.03 (82.1 ) 97.58 (98.5f)

Atomic position in R2O2S using lattice vector unit are (1/3, 2/3, u) for Y (La) with u  0.28, (1/3, 2/3, v) for O with Cal. data evaluated by ABINIT code with LDA method by Mikami and Nakamura [8]. Exp. data from Ref. [10]. Data obtained by chemical reaction: Y2O3 (s) + H2S (s) ? Y2O2S (s) + H2O (s) in present work. Cal. band gap by Mikami and Oshiyama [6]. Exp. band gap from Ref. [6]. Exp. data reported by Morosin [26]. Data obtained by chemical reaction: La2O3 (s)+H2S (s) ? La2O2S (s) + H2O (s) in present work. Data obtained by chemical reaction: La2O3 (s) + H2S (g) ? La2O2S (s) + H2O (g) from Ref. [11]. Cal. band gap by Vali [9]. Exp. data reported by Struck and Fonger [27].

5.15 5.82

Band gap c

1.14 1.61g (1.14h)

4.2 (2.61d, 4.8e) 4.2 (2.7a, 2.91i, 4.8j)

v  0.63, and (0, 0, 0) for S, respectively.

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

157

Fig. 2. Chemical bonding motifs of R2O2S (R = Y, La) crystal strcuture. (a) Capped trigonal anti-prism of RO4S3 (C3v); (b) Distorted tetrahedral of OR4 (C3v); (c) Trigonal antiprism of SR6 (C3v). The blue, pink, and yellow balls represent R, O, S atoms, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

they are very similar to each other in terms of the order of energy. Using La2O2S as an example, we have calculated the eigenergies of O (2s, 2p), S (3s, 3p) and La (5s, 5p, 5d and 6s) atoms by solving the non-relativistic Schrödinger equation in Opium code [29]. In our calculations, many-body effects were approximated by PBE exchange-correlation functional [30]. For the atomic orbital energy from the lowest to highest with respect to vacuum level (0 eV), we find that 35.90 eV (La, 5s), 23.80 eV (O, 2s), 22.30 eV (La, 5p), 17.14 eV (S, 3s), 8.97 eV (O, 2p), 6.94 eV (S, 3p), 3.67 eV (La, 5d) and 3.40 eV (La, 6s). This compares very well with the positions of different bands shown in Fig. 3b for La2O2S. The charge density of La2O2S is analyzed on (1 0 0) crystallographic plane from three different aspects. For Y2O2S, similar results are not shown here for bravity. The following parameters are calculated from planar electron density. 1. Local Wigner-Seitz radius (rs):

 rs ¼

1=3 3 4pnðx; yÞ

ð4Þ

2. Reduced Electron Density Gradient (s):

s¼

jrnðx; yÞj

ð3p2 Þ

1=3 4=3 n

ð5Þ

3. Electron Localization Function (ELF):

ELF ¼

1 1 þ a2

ð6Þ

where a is given by

a¼

Fig. 3. Calculated band structure, total and angular-momentum projected density of states of Y2O2S and La2O2S using sX-LDA exchange-correlation functional. The black dotted lines represent the Fermi level.

atom. The dominant bonding mechanism between RAO or RAS atomic pairs are ionic. Covalent interaction has a minor effect in both Y2O2S and La2O2S structures. The band gap of R2O2S structure is charge transfer type, i.e., the conduction band is mainly attributed to unoccupied (n  1)d and ns states of R, and valence bands on the other hand are predominated by s + p states of S and O. This conclusion is also confirmed by computing the atomic orbital energies of isolated R (Y or La), O and S atoms, and comparing them to the relative positions of different bands in solid, it is found that

s  sw ½n sunif ½n

ð7Þ

In Eq. (7), the kinetic energies s, sW and sunif are given in Refs. [31,32]. Sun et al. [32] has used different a values to distinguish chemical bonds in the construction of exchange-correlation functional. Due to Eq. (6), we may also use ELF to identify the types of chemical bonds in crystal structures. The reduced density gradient (s) serves as a additional criteria to indicate the local uniformity of electron density. Meanwhile, the local Wigner-Seitz radius characterizes the accumulation of electron density in space directly. Our results are shown in Fig. 4 on (1 0 0) plane for La2O2S. Fig. 4a shows that most valence electrons are localized around the cores of pseudo-atoms. The electron density at O positions is larger than that of S and La atoms. This is mainly due to the smaller atomic radius and larger electronegravity value of O than La and S elements. In the interestitial regions of a crystal structure, the local Wigner-Seitz radius is large due the low electron density. The slightly decreasing of rs between La and O or S clearly indicates the sharing of electron density due to the polar covalent bonds.

158

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

Fig. 4. Properties of electron density on (1 0 0) crystallographic plane of La2O2S. (a) Wigner-Seitz radius (rs); (b) reduced electron density gradient (s); and (c) electron localization function (ELF). The positions of different atoms on the plane are indicated by their atomic symbols. The core region of S atom shows anomalies in the calculated rs and s because of the pseudopotential approximation.

From Fig. 4b, we conclude that the valence electron density on (1 0 0) plane is high non-uniform in La2O2S structure. The pseudo-cores of La, O and S atoms are always closely surrounded by large s regions. The large s value is usually consistent with the exponential decreasing of electron density at the tail of valence orbitals [33]. The interestitial region also has small s value, because both electron density n and its gradient »n are small. Similar to Wigner-Seitz radius, s is relatively small in the bonding direction of LaAO or LaAS atomic pair, surrounding the pseudo-cores of La, O and S atoms. It should be noted that the significant decreasing of s in the middle between La and O or S is due to symmetry

constraint of charge density. The obtained s distribution is typical for polar covalent bonds [34,16]. The electron localization function depicted in Fig. 4c distinguishes the interestitial, bonding and ionic core regions in a better way than rs and s. As implied in Refs. [32,34], ELF is vey small ( 0) in the interestitial region of a crystal structure. Its value is roughly unity for completely filled spherical orbital. Uniform electron gas like behavior is defined by ELF = 0.5. From Fig. 4c, we can observe that OAO and SAO bonds do not exist in La2O2S structure. Meanwhile, LaAS and LaAO bonds are clearly indicated by the large ELF regions between two atoms. At O and S sites, ELF is close to unity, indicating the formation of closed

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

159

valence shells in anions. In summary, R2O2S structure is dominated by strong polar LaAO and LaAS bonds, and the ionic bond is the main bonding mechanism. In order to study the pressure-induced structural changes, both the lattice parameters and atomic positions were fully relaxed under isotropic hydrostatic pressure from 0 to 36 GPa. The ratios a/a0 and c/c0 were calculated and shown in Fig. 5a, where a0 and c0 represent the zero temperature/pressure values. All ratios decrease continuously with the increasing of pressure. The lattice along [0 0 1] direction is more compressible than that along [1 0 0] or [0 1 0] direction, since the ratio c/c0 changes more rapidly with pressure than a/a0. At the highest pressure (36 GPa) considered in this work, c-axis length is shortened by about 7% and 9% for Y2O2S and La2O2S respectively, showing a strong dependence on pressure. This anisotropic lattice compressibility is related to the layered crystal structure where chemical bonds between [RAO] layer and S atoms are relatively weak. Meanwhile, Y2O2S exhibits smaller compressibility than La2O2S under compression. The normalized bond lengths with the increasing of pressure are shown in Fig. 5b. All bond lengths decrease under compression in Y2O2S and La2O2S structures. YAO and LaAO bonds aligned in [0 0 1] direction are more resistant to compression than all other chemical bonds. YAS and LaAS have the longest bond lengths in the corresponding crystal structure. As a result, they are softer than YAO and LaAO bonds. The response of different chemical bonds to the hydrostatic pressure indicates that strongest chemical bond in R2O2S (R = Y or La) structure is the RAO bond parallel to c-axis. 3.2. Anisotropic elasticity The elastic constants determine the response of the crystal to external pressure which play an important role for engineering application as structural materials. In this work, the second order elastic constants Cij (i, j = 1, 2, 3, 4, 5, 6) were calculated by stress versus strain method using following strain modes [35].

½eT ¼ g½ 1 0 0 0 0 0 

ð8Þ

½eT ¼ g½ 0 0 1 1 0 0 

ð9Þ

where ½eT represents the transpose of the strain matrix with Voigt notation, and g the magnitude of strain. The maximum amplitude for each strain pattern was set as 0.003. The results are listed in Table 2. Both oxysulfides have the largest values of C11 and C33 indicating they are stiff against uniaxial stains e1 and e3. We note that C11 is bigger than C33 implying the structures are stiffer on basal plane (x-y plane) than c-axis. This property is in agreement with previous discussions about anisotropy in lattice compressibility and chemical bonds. All independent elastic constants of Y2O2S are larger than those of La2O2S, indicating a stronger chemical bonding in former structure than latter one. The electronegativity values of Y, La and O on the Pauling scaling are 1.22, 1.10 and 3.44. It is anticipated that YAO bond is less polarized than LaAO bond. Therefore, the larger mechanical moduli and elastic constants are obtained for Y2O2S. In Table 2, the bulk and shear moduli of oxysulfides were calculated by the Voigt-Ruess-Hill approximation [36]. Then, the Young’s modulus and Poisson’s ratio can be estimated by [37]:

E¼

9BVRH GVRH 3BVRH þ GVRH

ð10Þ

r¼

3BVRH  2GVRH 2ð3BVRH þ GVRH Þ

ð11Þ

Fig. 5. Variation of relative lattice constant a/a0, c/c0 (a), and normalized bond lengths as a function of pressure (b) from 0 to 36 GPa for R2O2S (R = Y, La).

The bulk moduli of Y2O2S and La2O2S are 135.6 and 113.4 GPa, respectively, which are lower than other common second phases in steel and iron materials like Fe3C (226.8 GPa [38]), TiC (242 GPa [39]), and TiB2 (251.4 GPa [40]), and even lower than BCC-iron (174 GPa [41]). Then, we want to analyze the mechanical moduli in three dimensions as a function of crystallographic direction. The detailed procedure is described in Appendix A. Fig. 6a and b show the dependence of bulk modulus on different crystalline directions. The anisotropy in bulk modulus is clear, because the surface contours of them show deviation from perfect spherical shape. More details can be found from the projections of the bulk modulus on the (0 0 1), (1 0 0), and (1 1 0) planes. The bulk modulus on (0 0 1) plane is isotropy (Fig. 6c), while weak anisotropy on (1 0 0) or (1 1 0) planes is revealed from Fig. 6d and e; i.e., inside (1 0 0) or (1 1 0) planes, B[0 0 1]
AB ¼

BV  BR BV þ BR

ð12Þ

160

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

Table 2 Calculated second order elastic constants (Cij in GPa), Bulk modulus (B, in GPa), shear modulus (G, in GPa), Young’s modulus (E, in GPa), Poisson’s ratio (r) and theoretical hardness (HV in GPa) of R2O2S (R = Y, La). Compounds

a b c

Cij

B

B/ G

E

r

HV

75.6 (79.9a)

1.8

191.26 (197.7a)

0.26 (0.24a)

9.7

59.8

1.9

152.6

0.28

7.4 (7.4c)

G

C11

C12

C13

C33

C44

C66

BV

BR

BVRH

GV

GR

GVRH

Y2O2S

230.7 (222.5a)

100.8 (79.2a)

84.4 (82.2a)

221.4 (201.1a)

89.9 (80.4a)

65.0 (71.6a)

135.8 (125.9a)

135.5 (125.0a)

76.4 (79.9a)

74.7 (81.2a)

La2O2S

194.5

80.1

75.5

171.3

67.8

57.2

113.6

113.1

135.6 (125.5a, 142b) 113.4

60.5

59.1

Cal. data obtained by Li and Ahuja [5]. Cal. data obtained by Mikami and Oshiyama [6]. The Knoop microhardness (HK) of 780 kg/mm2 was provided by Tsay and Wang [42], and this value is close to 7.4 GPa when converting it into Vicker’s hardness (HV).

Fig. 6. The surface construction of the bulk modulus of Y2O2S (a) and La2O2S (b). The (0 0 1), (1 0 0) and (1 1 0) planar projections of the bulk modulus are also shown in (c), (d) and (e), respectively.

Table 3 Calculated Bulk modulus (B, in GPa) and Young’s modulus (E, in GPa) in [1 0 0], [0 0 1], [0 1 1] and ½1 1 1 directions of R2O2S (R = Y, La). Compounds

Y2O2S La2O2S

AG ¼

GV  GR GV þ GR

B

E

B[1 0 0]

B[0 0 1]

E[1 0 0]

E[0 0 1]

E[0 1 1]

E½1 1 1

431.77 372.01

363.53 288.45

174.62 146.24

178.34 129.75

210.01 171.64

208.15 167.48

ð13Þ

On the other hand, Ravindran et al. [44] provided some shearing anisotropic factors (A1, A2 and A3), which are defined as:

A1 ¼

4C 44 C 11 þ C 33  2C 13

ð14Þ

A2 ¼

4C 55 C 22 þ C 33  2C 23

ð15Þ

A3 ¼

4C 66 C 11 þ C 22  2C 12

ð16Þ

In Table 4, we observe AG are lightly larger than AB, which agree very well with the previous discussion. The indexes A1, A2 and A3 defined in Ref. [44] are based on the strong directional dependent shear modulus. In R2O2S structure, we have A1 = A2 – A3 since C11 = C22, C13 = C23 and C44 = C55. A1 and A2 represent the anisotropy in shear modulus on x-z and y-z crystallographic planes. By the definition, A3 is unity due to the Cauchy condition C66 = 1/2 (C11  C12). Thus, A3 is a trivial index for R2O2S structure.

161

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

Fig. 7. The surface construction of the Young’s modulus of Y2O2S (a) and La2O2S (b). The (0 0 1), (1 0 0) and (1 1 0) planar projections of the Young’s modulus are also shown in (c), (d) and (e), respectively.

Table 4 Calculated longitudinal, transverse and average sound velocities (in km/s), the Debye temperature (in K), universal anisotropic index (AU), percent anisotropies (AB and AG), and shear anisotropic factors (A1, A2 and A3) of R2O2S (R = Y, La). Compounds

½1 0 0v l

½0 1 0v t1

½0 0 1v t2

½0 0 1v l

½1 0 0v t1

½0 1 0v t2

vl

vt

vm

bD

AU

AB

AG

A1

A2

A3

Y2O2S

6.69 (7.06a) 5.78 (6.03a)

3.55 (3.51a) 3.13 (2.67a)

4.18 (4.21a) 3.41 (3.22a)

6.56 (5.99a) 5.43 (4.96a)

4.18 (4.24a) 3.41 (3.33a)

4.18 (4.24a) 3.41 (3.33a)

6.787 (7.0b) 5.77

3.838 (4.1b) 3.211

4.269 (4.5b) 3.576

508 (527b) 395

0.116

0.0011

0.0113

1.269

1.269

1

0.1229

0.0022

0.0117

1.263

1.263

1

La2O2S a b

Data obtained by linear fitting acoustic phonon branches in phonon spectra performed by present work. From Li and Ahuja [5].

Ranganathan and Ostoja-Starzewski [45] proposed another universal anisotropy index (AU) as follows:

AU ¼ 5

GV BV þ 6P0 GR BR

ð17Þ

The results are shown in Table 4. A non-zero value implies elastic anisotropy. Obviously, the elastic anisotropy of La2O2S is stronger than that of Y2O2S. In addition, the ratio of B/G is widely used to reveal the ductility of materials. It is supposed that the B/G value of brittle material is lower than 1.75, and for ductile structure, the value is higher than 1.75 [46]. The obtained B/G values clearly imply that R2O2S (R = Y, La) are ductile phases. Experimentally, R2O2S could be formed by adding rare-earth elements in the liquid steel in order to reduce the oxide and sulfide impurities; and the fact is inspiring that the R2O2S products may not deteriorate the mechanical properties of alloy steels because of their ductile nature. In Ref. [47], Yang et al. also claimed that the rare earth elements could modify the sulfides into the more corrosion resistant oxysulfides of rare earth, and decrease the inclusions fraction. Besides, these oxysulfides would prefer to stay at the grain boundaries, which could prevent

the enrichment of MnS and FeS, and make the boundaries clean [48]. The hardness of second phase also plays an important role in steel making industry. Many theoretical models for estimating hardness of materials have been proposed. In this work, we employed a popular semi-empirical model proposed by Chen et al. [49]: Hv = 2(k2G)0.585  3, where k = G/B. The results are shown in Table 2. The theoretical hardness of Y2O2S is higher than La2O2S. The hardness of La2O2S is calculated to be 7.4 GPa, which is the same magnitude as the experimental data reported by Tsay and Wang [42]. 3.3. Anisotropy of acoustic velocities The anisotropy in sound velocity is obtained by solving the Christoffel equation for trigonal symmetry. Specifically, we have used Eq. (18) in the calculation [50].

  C11  qv 2    C21   C31

C12

C22  qv 2 C32

    ¼0  2 C33  qv

C13 C23

ð18Þ

where Cij = Cji, and all independent Cij components are written as

162

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

8 C11 ¼ C 11 n21 þ C 66 n22 þ C 44 n23 > > > > > C12 ¼ ðC 12 þ C 66 Þn1 n2 > > > < C13 ¼ ðC 13 þ C 44 Þn1 n3 2 2 2 > > C22 ¼ C 66 n1 þ C 11 n2 þ C 44 n3 > > > > > C23 ¼ ðC 13 þ C 44 Þn2 n3 > : C33 ¼ C 44 ðn21 þ n22 Þ þ C 33 n23

pffiffiffiffiffiffiffiffiffiffiffiffiffi ½1 0 0 : ½1 0 0v l ¼ C 11 =q; pffiffiffiffiffiffiffiffiffiffiffiffiffi ½0 0 1v t2 ¼ C 44 =q ð19Þ

here, Cij denotes elastic constants, and the directional cosines are given by n1 = sin h cos u, n2 = sin h sin u, and n3 = cos u in spherical coordinates. Besides some high-symmetry crystallographic directions ([1 0 0], [1 1 0], [0 0 1] and etc.), the eigenvalues and eigenvectors of Eq. (18) in an arbitrary direction are complex numbers. The complex solutions imply the attenuation of sound waves inside the bulk materials. In addition, three eigenvalues are obtained from Eq. (18), however, they do not represent the pure transverse and longitudinal wave modes. Generally, the one with largest eigenvalue is referred to faster mode. The remaining two eigenvalues are known as slower modes. In Fig. 8, the three-dimensional (3D) contour plots of three sound modes (faster mode (v þ 1 ) and  two slower modes (v  2 and v 3 ) are illustrated for Y2O2S structure. Similar plots are omitted here for La2O2S. All three sound modes are isotropic at the basal plane (x-y plane), and the weak anisotropy is seen for all modes on (1 0 0) crystallographic plane. The results are consistent with the crystallographic symmetry and also the calculated elastic constants of two structures. Although, the crystal structure of R2O2S (R = Y, La) consists of [R2O2]2+ and S2 alternative layers in [0 0 1] direction, the inter-layer bonds between R and S atoms are not van der Waals type weak interactions. From the obtained uniaxial elastic constants C11 and C33, we can see for both structures, the difference between them is small. The observed anisotropy in elasticity and sound velocity might be mainly caused by the additional RAO bond aligned in [0 0 1] direction. The band structure of sound velocity is displayed in Fig. 9 for (1 0 0) plane. The angle between [0 0 1] direction and another arbitrary crystallographic direction on (1 0 0) plane is defined by h. The large band gap is clearly observed between faster and other two slower modes. Similar to the band gap of electronic structure, the propagation of sound waves are forbidden in the band gap. Another interesting feature indicated in Fig. 9 is the faster mode and two slower modes are always out-of phase. Meanwhile, the latter two sound waves are in-phase. All sound modes in Y2O2S propagate faster than those of La2O2S due to the large density of latter compound. Solving Eq. (18) in high-symmetry crystallographic directions, the pure longitudinal and transverse modes are obtained [51,52]:

pffiffiffiffiffiffiffiffiffiffiffiffiffi ½0 0 1 : ½0 0 1v l ¼ C 33 =q; pffiffiffiffiffiffiffiffiffiffiffiffiffi ½0 1 0v t2 ¼ C 44 =q

½0 1 0v t1 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðC 11  C 12 Þ=2q; ð20Þ

½1 0 0v t1 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffi C 44 =q; ð21Þ

Here, vl represents the longitudinal sound velocity, vt1 the first transverse mode and vt2 the second transverse mode. The sound velocity anisotropy reveals the elastic anisotropy in the cell. The average (vm), the transverse (vl) and longitudinal (vt) sound velocity can be calculated from elastic moduli [53,54]:

vm ¼ vl ¼

   1=3 1 2 1 þ 3 v 3t v 3l

s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 4 1 ; Bþ G 3 q

vt ¼

ð22Þ pffiffiffiffiffiffiffiffiffi G=q

ð23Þ

Based on the computed longitudinal and transverse acoustic velocities, the Debye temperature bD is obtained by [55]:

HD ¼

   1=3 h 3n NA q vm k B 4p M

ð24Þ

Here h, kB, and NA are the Planck, Boltzmann and Avogadro constants, respectively. n refers to the number of atoms per formula unit, M is the mean molecular weight, and q is the density. The results are tabulated in Table 4. The obtained wave velocities are similar to other common second phases like Fe3C (around 3.7– 6.7 km/s [56]) and W2C (3.2–5.9 km/s [57]) in alloy steels. Meanwhile, it is known that bD can be used to characterize the strength of chemical bonds in the solids. The Debye temperature of Y2O2S is higher than La2O2S, so the chemical bonds among atoms in former structure are stronger than those in latter one. The results are consistent with the previous discussions about the electronic structures of them. 3.4. Anisotropy of thermodynamic properties The phonon dispersion and corresponding phonon DOS provide important information about the lattice dynamical properties of crystal. Many thermal-physical properties of solids depend on their phonon properties, e.g., specific heat, thermal expansion coefficient, heat conductivity. The calculation results by compress and expansion (±10%) of equilibrium crystal structure are shown in Fig. 10, and the phonon spectra for equilibrium structures can be found in Fig. 11. No imaginary frequencies can be found in the

  Fig. 8. Three-dimensional contour plots of three sound wave modes of Y2O2S. (a) Faster-mode (v þ 1 ); (b) slower-mode (v 2 ); and (c) slower-mode (v 3 ).

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

163

Fig. 9. Band structures of sound wave modes on (1 0 0) crystallographic plane. (a) Y2O2S and (b) La2O2S.

whole Brillouin zone for both oxysulfides proving their structure stability under certain lattice deformations. The crystal structure of R2O2S contains 5 atoms, giving rise to a total number of 15 phonon branches. At each k point, there will be 3 acoustic modes and 12 optical modes. As can be seen from Fig. 10, due to the significant mass contrast of R, O and S atoms in R2O2S structures, each sublattice vibrates almost independently. For example, the acoustic and some low-frequency optical phonons are mainly attributed to heavy R element, and those high-frequency optical modes are dominated by O-sublattice. Meanwhile, phonon branches involving S-sublattice are situated between phonon bands of R and O atoms. Compared with compressed structure, the phonon dispersions shift to lower frequency under the expansion of cell volume, since the force constants between atomic pairs in the crystal structure are usually reduced. Furthermore, acoustic branches in C-M ([0 1 0]) and C-A ([0 0 1]) directions of the both oxysulfides can be employed to estimate the longitudinal and transverse sound velocities directly [52]. The results are shown in Table 4, and the computed sound velocities are in good agreement with the values obtained by elastic constants. Meanwhile, we can also estimate elastic constants from sound velocities by Eqs. (20) and (21). The calculated elastic constants of Y2O2S are C11 = 256 GPa, C12 = 129 GPa, C33 = 184 GPa, C44 = 92 GPa; and the values of La2O2S are C11 = 211 GPa, C12 = 128 GPa, C33 = 143 GPa, C44 = 65 GPa. These results are comparable with the elastic constants calculated by stress-strain method (Table 4). The Grüneisen parameters of a crystal structure indicate its anharmonic effects in the phonon spectrum due to the change of cell volume [58]. In this paper, the phonon mode-Grüneisen constant (c (xi)) and macroscopic Grüneisen parameter (c) were calculated. The former one describing the shifting of phonon frequency with respect to the volume can be computed by the following expression:

cðxi Þ ¼ 

V 0 dxi ðVÞ dV

x0;i

ð25Þ

where xi represents the phonon mode i with angular frequency x, and V the cell volume. x0,i and V0 refer to the values of equilibrium structure. The macroscopic Grüneisen parameter is computed by:

P

cðTÞ ¼

cðx ÞC ðx ; TÞ P i V i i C V ðxi ; TÞ

i

here, the mode-specific heat is calculated by

ð26Þ

C V ðxi ; TÞ ¼ kB

 
hxi 2 expð
hxi =kB TÞ kB T ½expð
hxi =kB TÞ  12

ð27Þ

By expanding or compressing the equilibrium volume by 1%, the mode-Grüneisen parameters were calculated. We plot the mode-Grüneisen parameters shown in Fig. 11. Most c (xi) values are positive throughout the Brillouin zone for all phonon branches. The largest positive mode-Grüneisen parameter is given by a k-point in C-K direction. The eigenvector of this particular phonon mode is depicted in Fig. 12a. We can see that all atoms roughly vibrate in-phase in [0 0 1] direction. Moreover, La (Y) atoms show a larger displacement than those of O and S atoms in the same direction. The main reason why its frequency is very sensitive to the change of cell volume might be attributed to the stretching of La (Y)AO bond aligned in c-axis. Meanwhile, we also notice few negative values in the acoustic branches around L point. The phonon eigenvector corresponding to the one with the most negative value is illustrated in Fig. 12b. All atomic displacements are aligned in [1 0 0] direction on x-y basal plane. The phonon eigenvector represents an out-of-phase movement between [R2O2]2+ and S2 layers. The observed negative c (xi) seems to be relevant to the long-pair of chalcogen anions. Because the electrostatic interaction between OAS atoms is mainly repulsive, and the reduction is expected as a result of volume expansion. The overall macroscopic Grüneisen parameters of R2O2S structures is positive. Therefore, the thermal expansion of them shows normal behavior, i.e., volumetric and linear TECs increase with the increasing of temperature. Our results shown in Fig. 11 also indicate that the acoustic and some low-frequency optical phonons, which are mainly related to the motion of R (Y or La) and S atoms, are sensitive to the change of cell volume. Thus, they have large mode-Grüneisen parameter. In contrast, mode-Grüneisen parameter is small for high-frequency optical modes involving Osublattice only. The macroscopic Grüneisen parameters are shown in Fig. 13. All Grüneisen parameters decrease continuously with the temperature due to the increasing of total specific heat, and also because the high-frequency optical modes have smaller mode-Grüneisen parameters than acoustic modes (Fig. 11). For both oxysulfides, the Grüneisen parameter of c-axis is higher than that of basal plane (x-y plane). As a result, we would expect the linear TEC perpendicular to 2-D [RAO] layer will be higher than that in [1 0 0] direction. The thermal expansion coefficient, which characterizes the anharmonicity of cell structure, are very important for engineering

164

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

Fig. 10. The phonon dispersions and phonon density of states for Y2O2S (a) and La2O2S (b) at different volumes. The solid line is for the structure at compressed state (10%), while dot line represents the structure at expanded state (+10%) compared with the equilibrium volume.

materials. Based on the volumetric TEC (k) and macroscopic Grüneisen parameters, the linear TECs (a) can be calculated [58]:

kðTÞ ¼ 2aa ðTÞ þ ac ðTÞ

ð28Þ

aa ðTÞ c½1 0 0 ðTÞ ¼ ac ðTÞ c½0 0 1 ðTÞ

ð29Þ

where the ratio is computed from the flat region in Fig. 13 at high temperature. The calculated volumetric and linear TECs for R2O2S are shown in Fig. 14. The TECs for both oxysulfides are similar in temperature range 0–1500 K. At room temperature (300 K), the volumetric TECs for Y2O2S and La2O2S are 2.9
105 K1 and 2.8
105 K1, respectively. Usually, the thermal expansion mismatch between second phases and alloy steel matrix should be reduced as small as possible in order to improve their thermal shock resistance property, and thereby delaying the initiation and propagation of micro-crack at their interface. We note that for alloy steels, the average linear TEC is around 1
105 K1. The linear TEC of R2O2S has similar values, which may guarantee the thermal stability of alloy steels serving at different temperatures. Unfortunately, to our knowledge, the experimental TEC for R2O2S is not available at the moment. We still would like to compare our results to

experimental TEC data of another rare earth oxide, i.e. La2O3 with same space group. Wei et al. [59] provided experimental linear TECs of La2O3 between 300 and 2000 K. As can be seen from Fig. 14, the linear TECs of R2O2S are fairly close to the average linear TEC of La2O3. The dependence of thermal conductivity on temperature was calculated by the Slack’s method [60] between 300 and 1500 K. The behaviors of thermal transportation in this temperature range are considered to be relevant to their heat treatment process and hot rolling process.

k¼A

M H3D d n2=3 c2 ðTÞT

ð30Þ

where A is a constant, which is calculated from macroscopic Grüneisen parameters: A = 2.43
106/(1  0.514/c + 0.228/c2), M is the average mass per atom in the crystal, d is the cube root of the average volume per atom, n is the number of atoms in the primitive unit cell, and c (T) represents the macroscopic Grüneisen parameters. The calculated results are shown in Fig. 15. We can see clearly thermal conductivity decreases with the increasing of temperature by the law 1/T. La2O2S has lower thermal conductivity than Y2O2S for the obvious reasons, i.e., the heavier atomic mass of La than that

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

165

Fig. 11. The phonon dispersions and phonon mode-Grüneisen parameters for Y2O2S (a) and La2O2S (b) along high-symmetry directions at equilibrium volume. The color map is consistent with the increasing of phonon frequency. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Two special phonon eigenvectors are displayed for La2O2S phase. (a) Phonon mode with the largest mode-Grüneisen parameter (x = 2.22 THz, c = 5.28) at k-point (0.03, 0.06, 0.00) in C-K direction; (b) phonon mode with the most negative mode-Grüneisen parameter (x = 16.02 THz, c = 0.37) at L point (0.0, 0.5, 0.5). The module of eigenvector is indicated by the length of arrow. VESTA program is used for the visualization of eigenvectors.

166

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

Fig. 13. The macroscopic Grüneisen parameters for Y2O2S (a, b) and La2O2S (c, d) in two principal directions. (a, c) [1 0 0] direction; (b, d) [0 0 1] direction.

Fig. 14. Thermal expansion coefficients (TECs) as a function of temperature for R2O2S (R = Y, La). Linear TECs along the [0 0 1] direction are larger than that along [1 0 0] direction.

of Y leads to smaller phonon frequency and sound velocity in former structure. Finally, we note that our calculated average thermal conductivity value of La2O2S (5.3 W/m K in the temperature range of 300–1500 K) agrees very well with the experiment value reported by Tsay and Wang (5 W/m K) [42].

4. Conclusions Using first-principles calculations based on DFT, the structural stabilities, bonding characters, anisotropic mechanical and thermodynamic properties of R2O2S (R = Y, La) were investigated. The calculated band gap of 4.2 eV for both Y2O2S and La2O2S are close to the experimental data (4.8 eV); the dominant bonding mechanism of both structures are ionic. Under hydrostatic pressure, the normalized bond length of RAS bond decreases rapidly because of the relatively weak chemical bond between [R2O2]2+ and S2 layers. The mechanical modulus of R2O2S was evaluated from the elas-

Fig. 15. Thermal conductivities (k) as a function of temperature for R2O2S (R = Y, La).

tic constants. It was found that the Young’s modulus of R2O2S shows stronger anisotropy than bulk modulus; the maxima of Young’s modulus are found at the [0 1 1] and ½1 1 1 directions. The universal anisotropy index indicates that the elastic anisotropy of La2O2S is higher than that of Y2O2S. The theoretical hardness of R2O2S is also estimated, which is the same magnitude as experimental value. The Grüneisen parameters were computed based on phonon spectra. For both structures, values of c-axis are larger than those of basal plane (x-y plane). The volumetric TECs were calculated, and at room temperature, the values of Y2O2S and La2O2S are 2.9
105 K1 and 2.8
105 K1, respectively. The average thermal conductivities of Y2O2S and La2O2S are 10.1 and 5.3 W/m K between 300 and 1500 K, respectively. Acknowledgments This work was supported by the National Natural Science Foundation of China (51501139), the Science and Technology Project of Guangdong Province in China (Special Foundation for Additive

Y.F. Li et al. / Computational Materials Science 125 (2016) 154–167

Manufacturing Technologies: No. 2015B010122003, and Special Foundation for Practical Science and Technology Research and Development Project: No. 2015B090926009), the Postdoctoral Science Foundation funded project of China (No. 2014M552434). Appendix A. Graphical representation of the crystallographic dependence of the elastic moduli The directional dependence of the bulk modulus and the shear modulus for trigonal crystals can be calculated using the following relationship [61]:

1 2 ¼ ðS11 þ S12 þ S13 Þ  ðS11 þ S12  S13  S33 Þl3 B 1 ¼ ð1  E

2 2 l3 Þ S11

4

2

2

ðA1Þ 2

2

þ l3 S33 þ l3 ð1  l3 Þð2S13 þ S44 Þ þ 2l2 l3 ð3l1  l2 ÞS14 ðA2Þ

where S11, S12, etc. are elastic compliance constants, l1, l2 and l3 are calculated by cos a, cos b, and cos c, respectively. The angles a, b, and c are set with respect to the [1 0 0], [0 1 0] and [0 0 1] directions. Substituting the relationships of the direction cosines in spherical coordinates with respect to h and u (cos a ¼ sin h cos u, cos b ¼ sin h sin u, cos c ¼ cos h) into Eqs. (A1) and (A2) we obtain the equations used to plot the anisotropic mechanical figures shown in Figs. 6a_b and 7a_b. For trigonal crystal class, the uniaxial bulk moduli in [1 0 0], [0 1 0], [0 0 1] and [1 1 0] directions are expressed as

B½1 0 0 ¼ B½0 1 0 ¼ B½1 1 0 ¼

B½0 0 1

1 S11 þ S12 þ S13

1 ¼ 2S13 þ S33

ðA3Þ ðA4Þ

Moreover, the expressions for Young’s modulus in [1 0 0], [0 1 0], [0 0 1], [1 1 0] and [1 1 1] directions are given as:

E½1 0 0 ¼ E½0 1 0 ¼ E½1 1 0 ¼

1 S11

ðA5Þ

[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]

E½0 0 1 ¼

1 S33

ðA6Þ

E½0 1 1 ¼

4 S11 þ S33 þ 2S13 þ S44  2S14

ðA7Þ

E½1 1 1 ¼

9 ð2S11 þ S33 Þ þ 2ðS12 þ 2S13 Þ þ ð2S44 þ S66 Þ  4S14

ðA8Þ

References [1] D. Cavouras, I. Kandarakis, G.S. Panayiotakis, E. Kanellopoulos, D. Triantis, C.D. Nomicos, Appl. Radiat. Isotopes 49 (1998) 931. [2] A.A. Kaminskii, Laser Crystals, Springer-Verlag, Berlin, 1981. [3] Z. Yu, Y. Zhu, Rare Earth in Steel, Metallurgical Industry Press, Beijing, 1982.

[45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61]

167

L. Yue, L. Wang, J. Han, J. Rare Earth 28 (2010) 952. S. Li, R. Ahuja, J. Appl. Phys. 97 (2005) 103711. M. Mikami, A. Oshiyama, Phys. Rev. B 57 (1998) 8939. M. Mikami, A. Oshiyama, J. Lumin. 87–89 (2000) 1206. M. Mikami, S. Nakamura, J. Alloys Compd. 408–412 (2006) 687. R. Vali, Comp. Mater. Sci. 37 (2006) 300. M. Itoh, Y. Inabe, Phys. Rev. B 68 (2003) 035107. J.M. Polfus, T. Norby, R. Bredesen, J. Phys. Chem. C 119 (2015) 23875. B. Deng, G.H. Chan, F.Q. Huang, D.L. Gray, D.E. Ellis, R.P. Van Duyne, J.A. Ibers, J. Solid State Chem. 180 (2007) 759. M.D. Segall, P.J.D. Lindan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J. Clark, M.C. Payne, J. Phys-Condens. Mat. 14 (2002) 2717. A.E. Mattsson, P.A. Schultz, M.P. Desjarlais, T.R. Mattsson, K. Leung, Model. Simul. Mater. Sc. 13 (2005) R1. L.Z. Cao, J. Shen, N.X. Chen, J. Alloy Compd. 336 (2002) 18. P. Haas, F. Tran, P. Blaha, K. Schwarz, R. Laskowski, Phys. Rev. B 80 (2009) 195109. R. Gillen, S.J. Clark, J. Robertson, Phys. Rev. B 87 (2013) 125116. J.L.F. da Silva, M.V. Ganduglia-Pirovano, J. Sauer, V. Bayer, G. Kresse, Phys. Rev. B 75 (2007) 045121. D.M. Bylander, L. Kleinman, Phys. Rev. B 41 (1990) 7868. A. Seidl, A. Gorling, P. Vogl, J.A. Majewski, M. Levy, Phys. Rev. B 53 (1996) 3764. S.J. Clark, J. Robertson, Phys. Rev. B 82 (2010) 085208. H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. S. Baroni, S. de Gironcoli, A.D. Corso, P. Giannozzi, Rev. Mod. Phys. 73 (2001) 515. C. Lee, X. Gonze, Phys. Rev. B 51 (1995) 8610. G. Kern, G. Kresse, J. Hafner, Phys. Rev. B 59 (1999) 8551. B. Morosin, Acta Cryst. B29 (1973) 2647. C.W. Struck, W.H. Fonger, Phys. Rev. B 4 (1971) 22. J.W. Haynes, J.J. Brown, J. Electrochem. Soc. 115 (1968) 1060. A.M. Rappe, K.M. Rabe, J.D. Joannopoulos, Phys. Rev. B 41 (1990) 1227. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. J. Tao, J.P. Perdew, V.N. Staroverov, G.E. Scuseria, Phys. Rev. Lett. 91 (2003) 146401. J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G.I. Csonka, G.E. Scuseria, J.P. Perdew, Phys. Rev. Lett. 111 (2013) 106401. J. Sun, B. Xiao, A. Ruzsinszky, J. Chem. Phys. 137 (2012) 051101. B. Xiao, J. Sun, A. Ruzsinszky, J. Feng, R. Haunschild, G.E. Scuseria, J.P. Perdew, Phys. Rev. B 88 (2013) 184103. J.J. Zhao, J.M. Winey, Y.M. Gupta, Phys. Rev. B 75 (2007) 094105. W. Zhou, L. Liu, B. Li, P. Wu, Q. Song, Comp. Mater. Sci. 46 (2009) 921. H. Xiang, Z. Wu, Inorg. Chem. 47 (2008) 2706. J.H. Jang, I.G. Kim, H.K.D.H. Bhadeshia, Comp. Mater. Sci. 44 (2009) 1319. J.J. Gilman, B.W. Roberts, J. Appl. Phys. 32 (1961) 1405. L. Sun, Y. Gao, B. Xiao, Y. Li, G. Wang, J. Alloys Comp. 579 (2013) 457. H.C. Herper, E. Hoffmann, P. Entel, Phys. Rev. B 60 (1999) 3839. B.J. Tsay, I.H. Wang, Mater. Lett. 34 (1998) 180. D.H. Chung, W.R. Buessem, Anisotropy in Single Crystal Refractory Compounds, Plenum Press, New York, 1968. P. Ravindran, L. Fast, P.A. Korzhavyi, B. Johansson, J. Wills, O. Eriksson, J. Appl. Phys. 84 (1998) 4891. S.I. Ranganathan, M. Ostoja-Starzewski, Phy. Rev. Lett. 101 (2008) 055504. S.F. Pough, Philos. Mag. 45 (1954) 823. H.J. Yang, K. Yang, B.C. Zhang, Mater. Lett. 61 (2007) 1154. L.M. Wang, Q. Lin, J.W. Ji, D.N. Lan, J. Alloys Comp. 408–412 (2006) 384. X.Q. Chen, H. Niu, D. Li, Y. Li, Intermetallics 19 (2011) 1275. J. Feng, B. Xiao, R. Zhou, W. Pan, Acta Mater. 61 (2013) 7364. K. Brugger, J. Appl. Phys. 36 (1965) 768. G. Grimval, Thermophysical Properties of Materials, North Holland, 1999. Y. Ding, B. Xiao, Comp. Mater. Sci. 82 (2014) 202. O.L. Anderson, J. Phys. Chem. Solids 24 (1963) 909. C. Kittel (Ed.), Introduction to Solid State Physics, third ed., John Wiley, 1966. W.H. Zhang, Z.Q. Lv, Z.P. Shi, S.H. Sun, Z.H. Wang, W.T. Fu, J. Magn. Magn. Mater. 324 (2012) 2271. Y.Z. Liu, Y.H. Jiang, R. Zhou, J. Feng, Ceram. Int. 40 (2014) 2891. Y.C. Ding, B. Xiao, RSC Adv. 5 (2015) 18391. C. Wei, J.L. Fan, H.R. Gong, J. Alloys Comp. 618 (2015) 615. G.A. Slack, Solid State Phys. 34 (1979) 1. J.F. Nye, Physical Properties of Crystals, Oxford University Press, Oxford, 1985.