I. INTRODUCTION
Ultra supercritical (USC), owing to high efficiency, low energy consumption and low emission, have been extensively used in thermal power industry. Reference Viswanathan and Bakker1–Reference Viswanathan, Sarver and Tanzosh3 The materials of 18-8 austenitic heat resistant steel, such as TP304H, TP347H, have been applied in above 600 °C of superheater/reheater in USC. Reference Li, Zou, Zhang, Zou and Du4,Reference Zhao, Xie and Smith5 However, the internal oxidation films generate large compressive stress in cooling, which result in oxidation films easily peel off at working. Reference Pantleon and Montgomery6 The oxidation films are Fe3O4 and (Fe, Cr)3O4. Moreover, the oxidation films may cause blockage, which leading to the pipeline burst. To date, some works have been reported as the ground state properties of Fe3O4 either by theoretical calculations or by empirical methods. The Curie temperature of bulk Fe3O4 is close to 860 K. Reference Masrour, Hlil, Hamedoun, Benyoussef, Mounkachi and Moussaoui7 The magnetite has been considered as a mixed-valence 3d transition-metal compound with the chemical formula FeA 3+[Fe2+Fe3+]BO4 (Ref. Reference Jeng, Guo and Huang8); the tetrahedral (A) sites are hold by Fe3+ and the octahedral (B) sites are occupied by equal numbers of Fe2+ and Fe3+, which result in an inverse spinel structure. Reference Wang, Lee, Zhang, Shang, Chen, Kovacs and Liu9,Reference Mohammadi, Barikani and Barmar10 Lin et al. Reference Lin, Chu and Wang11 have synthesized magnetite nanoparticles by mechanochemical reaction and then annealed at 200, 400, 600, and 800 °C for 1 h in 10% H2/Ar mixed gas. Results indicated that the saturation magnetization of Fe3O4 nanoparticles increased gradually with the increase of annealing temperature up to 800 °C. Iijima et al. Reference Sato, Iijima, Seki and Inagaki12 have synthesized ultrafine ferrite nanoparticles, such as Fe3O4, CoFe2O4, and MnFe2O4, by co-precipitation process, and also investigated the effect of annealing on the morphology and magnetic properties of these ferrites. Using first-principles methods, Odkhuu and co-workers Reference Odkhuu, Taivansaikhan, Yun and Hong13 have investigated the effects of charge ordering in the octahedral sites of Fe3O4 and CoFe2O4 on their magnetostrictions. Cheng et al. Reference Cheng, Li, Wang, Liu, Ren, Cui and Zheng14 have estimated the electrical and magnetic properties of Zn-doped Fe3O4 at different doping concentrations of Zn. The results showed that Zn x Fe3−x O4 (0 ≤ x ≤ 0.875) is half-metallic with a full spin polarization. A half-metal to metal transition at the (100) Fe3O4 surface has also been observed by means of spin-resolved photoemission experiments on epitaxial and high quality thin films as well as by density functional theory (DFT) with generalized gradient approximation (GGA). Reference Fonin, Pentcheva, Dedkov, Sperlich, Vyalikh, Scheffler, Rüdiger and Güntherodt15 Based on aforementioned discussions, the properties of Fe3O4 and (Fe, Cr)3O4 are rarely systematically investigated. Moreover, these properties are very important for the research of ultra supercritical. In this work, we use the precise first-principles to calculate the elastic anisotropy and anisotropy of thermal conductivity of Fe3−x Cr x O4 compounds avoiding the experiments. The results are crucial in understanding the instinct physical properties and providing a theoretical basis for obtaining excellent performance of Fe3−x Cr x O4 compounds.
II. METHODS AND DETAILS
First principles calculations based on density functional theory (DFT) were used to calculate the ground properties, as implemented in Cambridge Serial Total Energy Package (CASTEP) code. Reference Segall, Lindan, Probert, Pickard, Hasnip, Clark and Payne16 The ultrasoft pseudopotentials were used to evaluate the interactions between valence electrons and ionic cores. Reference Vanderbilt17 The valence electrons of atoms Cr, Fe, and O are considered as 3s 23p 63d 54s 1, 3s 23p 63d 64s 2, and 2s 22p 4, respectively. The exchange-correlation potential of geometry optimizations and property calculations were calculated using generalized gradient approximation (GGA) within the Perdew–Burke–Ernzerhof (PBE). Reference Perdew, Burke and Ernzerhof18 The equilibrium crystal structures of these compounds were obtained using BFGS (Broyden–Fletcher–Goldfarb–Shannon) optimized method. Reference Liu, Jiang, Feng and Zhou19 The electron energy minimization process was calculated using the Pulay density mixing scheme. The cut-off energy was 450 eV and the first irreducible Brillouin zone (BZ) integration was set as 8 × 8 × 8 Monkhorst–Pack Reference Monkhorst and Pack20 k-point sampling for the calculation of properties for the entire compounds. The total energy changes during the optimization process were finally converged to 1 × 10−8 eV, and the forces per atom were reduced to 0.01 eV/Ǻ.
The crystalline structure of Fe3O4 has the space group
$Fd\bar 3m$
(no. 227) with eight formula unit (Z = 8). The Wyckoff positions of atoms are site 8a FeI (0.125, 0.125, 0.125), site 16d FeII (0.5, 0.5, 0.5), and site 32e O (0.2539, 0.2539, 0.2539). The structures of Fe3−x
Cr
x
O4 (x = 0–3) compounds are shown in Fig. 1. When chromium atoms occupy the site 8a (FeI), it will get Fe2CrO4. If chromium atoms occupy the site 16d (FeII), it will get FeCr2O4. The stress–strain method was used to obtain elastic coefficient.
Reference Fast, Wills, Johansson and Eriksson21
In this method, the stress–strain curve was described by the Hooker's law:
where σ ij is the stress tensor; ε kl is the Lagrangian strain tensor; C ijkl is the elastic constant tensor which is a 6 × 6 matrix. Several different strain modes were imposed on the crystal structure, and the Cauchy stress tensor for each strain mode was evaluated. The crystal symmetry was used to estimate the elastic constant tensor; for a cubic crystal, the equation can be expressed as Reference Liu, Jiang, Xing, Zhou and Feng22 :
$$\left( {\matrix{ {{\sigma _1}} \cr {{\sigma _2}} \cr {{\sigma _3}} \cr {{\tau _4}} \cr {{\tau _5}} \cr {{\tau _6}} \cr} } \right) = \left( {\matrix{ {{C_{11}}} & {{C_{12}}} & {{C_{12}}} & 0 & 0 & 0 \cr {{C_{12}}} & {{C_{11}}} & {{C_{12}}} & 0 & 0 & 0 \cr {{C_{12}}} & {{C_{12}}} & {{C_{11}}} & 0 & 0 & 0 \cr 0 & 0 & 0 & {{C_{44}}} & 0 & 0 \cr 0 & 0 & 0 & 0 & {{C_{44}}} & 0 \cr 0 & 0 & 0 & 0 & 0 & {{C_{44}}} \cr} } \right)\left( {\matrix{ {{\varepsilon _1}} \cr {{\varepsilon _2}} \cr {{\varepsilon _3}} \cr {{\gamma _4}} \cr {{\gamma _5}} \cr {{\gamma _6}} \cr} } \right)\quad ,$$
where C ij is the elastic constants, σ i and τ i are normal stress and shearing stress, respectively, and ε i and γ i refer to corresponding normal strain and shearing strain, respectively. The required strain patterns are significantly reduced in the stress–strain method compared to the energy-strain method. Reference Ravindran, Fast, Korzhavyi, Johansson, Wills and Eriksson23,Reference Xiao, Feng, Zhou, Jiang and Zhou24 There are three different independent symmetry elements (C 11, C 12, and C 44) for cubic crystal.
The crystal structures of Fe3−x Cr x O4 (x = 0–3) compounds.
The mechanical properties of polycrystalline materials, such as bulk modulus (B), shear modulus (G), Young's modulus (E), and Poisson's ratio (σ), can be estimated using the Voigt–Reuss–Hill approximation. Reference Wu, Chong, Zhou, Jiang and Feng25 The formulas were as follows Reference He, Jiang, Zhou and Feng26 :
here, B V, B R, and B VRH are the bulk modulus calculated by the Voigt, Reuss, and Voigt–Reuss–Hill approximation method, respectively. G V, G R, and G VRH are the shear modulus estimated within the Voigt, Reuss, and Voigt–Reuss–Hill approximation method, respectively. In addition, the compressibility (λ) and shear stiffness (μ) can be also estimated using the following expressions:
III. RESULTS AND DISCUSSIONS
A. Mechanical properties
The mechanical stability of a crystal implies that the strain energy must be positive. The intrinsic mechanical stability criteria of a cubic crystal using Born–Huang's lattice dynamical theory as follows Reference Liu, Jiang, Feng and Zhou19 :
As displayed in Table I, the calculated elastic constants of Fe3−x Cr x O4 compounds satisfy the aforementioned criteria, illustrating that these compounds are mechanically stable structures. C 11 represents the uniaxial deformation along [001] direction. C 12 is a mixture of σ11 and σ22 or the pure shear stress at (110) crystal plane along [110] direction. Reference He, Jiang, Zhou and Feng26 In Fe3−x Cr x O4 compounds, Cr3O4 has the largest C 11 value with 308.3 GPa, indicating that Cr3O4 has the highest incompressibility under uniaxial stress along crystallographic a axis (ε11) among Fe3−x Cr x O4 compounds. C 44 refers to the shear modulus on (100) crystal plane. The C 44 values of Fe3O4, Fe2CrO4, FeCr2O4, and Cr3O4 are 35.6 GPa, 59.7 GPa, 36.1 GPa, and 56.3 GPa, respectively.
The calculated lattice parameters, elastic constants (C ij ), bulk modulus (B), shear modulus (G), Young's modulus (E) and Poisson's ratio (σ) of Fe3−x Cr x O4 (x = 0, 1, 2, 3) (in GPa).

a Cal. in Ref. Reference Hu, Gao and Yang27.
b Cal. in Ref. Reference Roldan, Carballal and Leeuw28.
c Exp. in Ref. Reference Reichmann and Jacobsen29.
Bulk modulus reflects the compressibility of a solid under hydrostatic pressure. According to the Voigt–Reuss–Hill approximation, the bulk modulus (B), shear modulus (G), Young's modulus (E), Poisson's ratio (σ), the first order Lame constant (λ) and the second Lame constant (μ) are summarized in Table I. Moreover, the calculated results are consistent with the previously available results. Reference Hu, Gao and Yang27–Reference Reichmann and Jacobsen29 The calculated bulk modulus, shear modulus, and Young's modulus of these compounds as a function of Cr doped atomic number are plotted in Fig. 2(a). As shown in Fig. 2(a), with the increase of doping chromium atomic number, the moduli of these compounds firstly decrease, and then increase. The hardness of a material may be more related to shear modulus based on the work of Jiang and co-workers. Reference Jiang, Zhao and Jiang30 The larger the shear modulus (G) is, the higher the hardness of a compound. Therefore, the sequence of the hardness of Fe3−x Cr x O4 compounds: Cr3O4 > Fe2CrO4 > Fe3O4 > FeCr2O4. The Poisson's ratio is an indicator for the degree of directionality of chemical bonds, which is σ = 0.1 for covalent bonds and σ = 0.25 for ionic bonds. Reference Korozlu, Colakoglu, Deligoz and Aydin31 The calculated Poisson's ratio of Fe3−x Cr x O4 compounds ranges from 0.34 (Fe2CrO4) to 0.42 (FeCr2O4), which clearly indicates complex bond in Fe3−x Cr x O4 compounds. One of the important parameters is the ductile/brittle behavior, which is described by both classical criteria of Pugh's modulus ratio B/G (Ref. Reference Pugh32) and Cauchy pressure C 12 − C 44. Reference Fu and Yoo33 A compound is consistent with a ductile property when B/G is larger than 1.75; otherwise, it is a brittle property. Moreover, the positive value (C 12 − C 44) implies a metallic bonding based on Pettifor's work. Reference Pettifor34 Figure 2(b) shows the calculated B/G and (C 12 − C 44) as a function of Cr doped atomic number. As presented in Fig. 2(b), it can be found that B/G varies from 2.86 (Fe2CrO4) to 5.83 (FeCr2O4), which implies that all of compounds are ductile.
(a) Calculated bulk modulus (B), Young's modulus (E) and shear modulus (G) of Fe3−x Cr x O4. (b) Pugh's modulus ratio B/G and Cauchy pressure C 12 − C 44 of Fe3−x Cr x O4.

It is needed to characterize the extent of anisotropy because all single crystals are anisotropic. However, many different ways can be used to evaluate the elastic anisotropy. The universal anisotropic index (A U) Reference Ranganathan and Starzewski35 and the percent anisotropy (A B and A G) Reference Ravindran, Fast, Korzhavyi, Johansson, Wills and Eriksson23 are used to characterize the elastic anisotropy in this work. Equations were as follows:
where B V, G V, B R, and G R are the bulk modulus and shear modulus evaluated using Voigt and Reuss methods, respectively. If all indexes are zero, it is an isotropic structure. The large discrepancies from zero illustrate highly mechanical anisotropic properties. The results are tabulated in Table II. The universal anisotropic index (A U) is an important factor to describe the degree of mechanical anisotropy of a compound. The larger the value of A U is, the stronger the anisotropy of a compound. Therefore, the elastic anisotropy order of Fe3−x Cr x O4 compounds: Fe2CrO4 > Fe3O4 > FeCr2O4 > Cr3O4 (Table II).
The calculated values for universal anisotropic (A U) and percent anisotropy (A B and A G).

Besides, the elastic anisotropy of crystalline structure is illustrated using the curved surface of three-dimensional (3-D) structure. Here, Young's modulus was plotted at different directions. The directional dependence of Young's modulus of cubic crystal can be calculated by Reference Liu, Jiang, Feng and Zhou19 :
$$E(\vec n) =\; {1 \over {{S_{11}} - {\beta _1}\left( {n_1^2n_2^2 + n_1^2n_3^2 + n_2^2n_3^2} \right)}},\quad {\rm{where\ }}{\beta _1} = 2{S_{11}} - 2{S_{12}} - {S_{44}}\quad ,$$
where S ij is the elastic compliance constant and n 1, n 2, and n 3 are the directional cosines. Substituting the relationships of the direction cosines in spherical coordinates with respect to θ and φ (n 1 = sin θ cos φ, n 2 = sin θ sin φ, n 3 = cos φ) into Eq. (13). The 3-D structures of Young's modulus of Fe3O4, Fe2CrO4, FeCr2O4, and Cr3O4 are shown in Fig. 3. In each graph, the surface configuration is the magnitude of Young's modulus on different orientations. From Fig. 3, there is a strong anisotropic of Young's modulus along different orientations, especially for Fe2CrO4 shows strong anisotropy because the surface contour deviated from a spherical shape. The projections of Young's modulus on different crystal planes are depicted in Fig. 4. The expressions of Young's modulus in the three principal directions of cubic structure are given as E [100] = E [010] = E [001] = 1/S 11, while in [100] direction, it is calculated by E [110] = 4/(2S 11 + S 12 + S 44). As observed in Fig. 4, at (001) plane, the graph of Cr3O4 is almost sphere, which result in almost isotropic. Fe3O4, Fe2CrO4, and FeCr2O4 show obviously anisotropic character. Moreover, the calculated Young's modulus of Fe2CrO4 is more anisotropic than other compounds at [110] directions.
The surface constructions of Young's modulus (E) for (a) Fe3O4, (b) Fe2CrO4, (c) FeCr2O4, and (d) Cr3O4, respectively. The units are in GPa.
The (a) (001) and (b) (110) plane projections of Young's modulus for Fe3−x Cr x O4 compounds.

B. Anisotropy of acoustic velocities
Debye temperature is related to many physical properties such as specific heat, elastic constant, and melting point. Reference Wang and Sun36 The Debye temperature (ΘD) can be calculated based on average sound velocity v m, the expressions were given by Reference Feng, Xiao, Chen, Du, Yu and Zhou37 :
$${v_{\rm{l}}} = \sqrt {{{\left( {B + \left( {{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-\nulldelimiterspace} 3}} \right)G} \right)} \over \rho }} \quad ,$$
$${v_{\rm{t}}} = \sqrt {{G \over \rho }}\quad ,$$
here, h is Planck's constant, k is Boltzmann constant, N A is Avagadro's constant, n is number of atoms in the molecule, M is molecular weight, and ρ is theoretical density of the compound. B and G refer to the bulk and shear modulus, respectively. v l and v t are the longitudinal sound velocity and transverse sound velocity, respectively. The calculated Debye temperatures and sound velocities of Fe3−x Cr x O4 compounds are listed in Table III. Debye temperature reflects the strength of chemical bonding in crystal structure, which inverses to the molecular weight. As observed in Table III, Cr3O4 has the largest chemical bonding among Fe3−x Cr x O4 compounds, because Cr3O4 has the largest Debye temperature with 507.6 K. In general, low density and high bulk modulus lead to relatively large sound velocity. Thus, Cr3O4 has the largest sound velocity among Fe3−x Cr x O4 compounds, which is excellent agreement with aforementioned discussion.
Sound velocities and Debye temperatures of Fe3−x Cr x O4 compounds.

The sound velocity of a solid also exist anisotropy, which is evaluated by the symmetry of crystal and propagation directions. Such as, the pure transverse and longitudinal modes can be only found for [001], [110], and [111] directions in cubic crystal, and the quasitransverse or quasilongitudinal waves are found in other directions. In different directions, the sound velocities of cubic crystal are correlated with elastic constants and phonon frequencies, which can be written as Reference Feng, Xiao, Chen, Du, Yu and Zhou37 :
$$\left[ {110} \right]{v_{\rm{l}}} = \sqrt {{{({C_{11}} + {C_{12}} + 2{C_{44}})} \mathord{\left/ {\vphantom {{({C_{11}} + {C_{12}} + 2{C_{44}})} {2\rho }}} \right. \kern-\nulldelimiterspace} {2\rho }}} ;\quad \left[ {1\bar 10} \right]{v_{{{\rm{t}}_1}}} = \sqrt {{{({C_{11}} - {C_{12}})} \mathord{\left/ {\vphantom {{({C_{11}} - {C_{12}})} \rho }} \right. \kern-\nulldelimiterspace} \rho }} ;\quad \left[ {001} \right]{v_{{{\rm{t}}_2}}} = \sqrt {{{{C_{44}}} \mathord{\left/ {\vphantom {{{C_{44}}} \rho }} \right. \kern-\nulldelimiterspace} \rho }}$$
$$\!\!\left[ {111} \right]{v_{\rm{l}}} = \sqrt {{{({C_{11}} + 2{C_{12}} + 4{C_{44}})} \mathord{\left/ {\vphantom {{({C_{11}} + 2{C_{12}} + 4{C_{44}})} {3\rho }}} \right. \kern-\nulldelimiterspace} {3\rho }}} ;\quad \left[ {11\bar 2} \right]{v_{{{\rm{t}}_1}}} = {v_{{{\rm{t}}_2}}} = \sqrt {{{({C_{11}} - {C_{12}} + {C_{44}})} \mathord{\left/ {\vphantom {{({C_{11}} - {C_{12}} + {C_{44}})} {3\rho }}} \right. \kern-\nulldelimiterspace} {3\rho }}}$$
where v
l is the longitudinal wave of sound velocity;
${v_{{{\rm{t}}_1}}}$
refers to the first transverse wave of sound velocity;
${v_{{{\rm{t}}_2}}}$
is the second transverse of sound velocity. The calculated sound velocities of these compounds are tabulated in Table IV. The anisotropic properties of sound velocities are related to the elastic anisotropy. For example, the C
11 determines the longitudinal sound velocities along [100] direction and C
44 corresponds to the transverse mode.
The anisotropic sound velocities of Fe3−x Cr x O4, and the unit of velocity (v) is m/s.

C. Minimum thermal conductivity and anisotropy
Thermal conductivity describes the diffusivity of heat flow by phonon transport in a temperature gradient. The value of minimum conductivity is important for high temperature applications, since the conductivity decreases with the increase of temperature to a limiting value known as the minimum thermal conductivity. Reference Feng, Xiao, Zhou, Pan and Clarke38 In this section, the minimum thermal conductivity was calculated by both the Cahill's model and Clarke's model, which can be estimated as Reference Clarke39,Reference Cahill, Watson and Pohl40 :
$${\rm{Clark&apos1;s\ model\!:\ }}{k_{\min }} = 0.87{k_{\rm{B}}}{\overline {{M_{\rm{a}}}} ^{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}{E^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}}{\rho ^{{1 \mathord{\left/ {\vphantom {1 6}} \right. \kern-\nulldelimiterspace} 6}}},\quad \overline {{M_{\rm{a}}}} = \left[ {{M \mathord{\left/ {\vphantom {M {(m \cdot {N_{\rm{A}}})}}} \right. \kern-\nulldelimiterspace} {(m \cdot {N_{\rm{A}}})}}} \right]\quad ,$$
here,
$\overline {{M_{\rm{a}}}}$
is average mass per atom, ρ is theoretical density of the compound, M refers to molar mass, m is total number of atoms per formula, N
A is Avogadro's constant, k
B is Boltzmann's constant, n is density of number of atoms per volume, v
l is the longitudinal sound velocity,
${v_{{{\rm{t}}_1}}}$
and
${v_{{{\rm{t}}_2}}}$
are the transverse sound velocity, and E is the Young's modulus. According to theoretical density ρ, Young's modulus E, longitudinal v
l and transverse sound velocity v
t, the minimum thermal conductivity could be deduced in Table V. From Table V, there is a little difference between the two models of minimum thermal conductivities, and the values obtained by the Cahill's model are slightly larger than that of by the Clarke's model. It is worth noting that the predicted minimum thermal conductivity is the intrinsic thermal conductivity of a perfect crystal. Actually, the phonon scattering and microstructure of a material will influence thermal conductivity. The directional dependence of minimum thermal conductivity can be depicted based on the anisotropy of Young's modulus. The 3-D surface contour of minimum thermal conductivity of Fe3−x
Cr
x
O4 compounds are plotted in Fig. 5. It is shown that all of Fe3−x
Cr
x
O4 compounds except for Cr3O4 have an anisotropic minimum thermal conductivity. Because the contour of Cr3O4 is almost sphere and it is verified by the anisotropy coefficient A
U = 0.02 (in Table II). Moreover, the anisotropic minimum thermal conductivity of a crystal at different crystal planes is also presented. For cubic crystal, the projections of minimum thermal conductivity at (001) and (110) crystal planes were given in spherical coordinates
Reference Feng, Xiao, Zhou and Pan41
:
$$\kappa _{\min }^{(001)} = 0.87{k_{\rm{B}}}{\overline {{M_{\rm{a}}}} ^{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}{\rho ^{{1 \mathord{\left/ {\vphantom {1 6}} \right. \kern-\nulldelimiterspace} 6}}}{\left( {{S_{11}}\left( {{{\sin }^4}\,\varphi + {{\cos }^4}\,\varphi } \right) + \left( {{S_{44}} + 2{S_{12}}} \right){{\cos }^2}\,\varphi \;{{\sin }^2}\,\varphi } \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}}$$
$$\kappa _{\min }^{(110)} = 0.87{k_{\rm{B}}}{\overline {{M_{\rm{a}}}} ^{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}{\rho ^{{1 \mathord{\left/ {\vphantom {1 6}} \right. \kern-\nulldelimiterspace} 6}}}{\left( {{S_{11}}\left( {{1 \over 2}{{\sin }^4}\,\theta + {{\cos }^4}\,\theta } \right) + \left( {{S_{44}} + 2{S_{12}}} \right)\left( {{{\cos }^2}\,\theta {{\sin }^2}\,\theta + {1 \over 4}{{\sin }^4}\theta } \right)} \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}}\quad .$$
The calculated thermal conductivity (k) of Fe3−x Cr x O4 compounds by Clarke's model and Cahill's model. The unit of the minimum thermal conductivity is W/(m K).

Contour plots of minimum thermal conductivity of Fe3−x Cr x O4 compounds in 3-D space.
The results are illustrated in Fig. 6. From the projection images, one can clearly see that minimum thermal conductivity of FeCr2O4 along [010] direction is smaller than that along other directions of other compounds and the peak value occurs along [110] direction at (001) plane, which are 1.36 W/(m K), 1.51 W/(m K), 1.24 W/(m K), and 1.63 W/(m K) for Fe3O4, Fe2CrO4, FeCr2O4, and Cr3O4, respectively. Moreover, minimum thermal conductivity of Cr3O4 along
$\left[ {1\bar 11} \right]$
direction is larger than that along the same direction of other compounds and minimum thermal conductivity of FeCr2O4 along
$\left[ {1\bar 10} \right]$
direction with 1.05 W/(m K) is the smallest at (110) plane.
Planar projections of minimum thermal conductivity of Fe3−x Cr x O4 compounds at (001) and (110) crystallographic planes.

IV. CONCLUSIONS
In summary, we have investigated the mechanical properties and thermal properties of Fe3−x Cr x O4 compounds using first-principles calculations. Elastic constants were evaluated using the stress–strain method. These compounds are mechanically stable structures, because they satisfy the mechanical stability criteria. The mechanical moduli were also evaluated and the calculated bulk moduli of Fe3O4, Fe2CrO4, FeCr2O4, and Cr3O4 are 190.9 GPa, 135.5 GPa, 180.1 GPa, and 235.6 GPa, respectively. The mechanical anisotropy was illustrated by both of anisotropic indexes and 3-D surface contour. The calculated Debye temperature of these compounds varies from 380.4 K (FeCr2O4) to 507.6 K (Cr3O4). Minimum thermal conductivity of Fe3−x Cr x O4 compounds was calculated by both the Clarke's model and Cahill's model. The calculated values of Fe3O4, Fe2CrO4, FeCr2O4, and Cr3O4 are 1.34 W/(m K), 1.06 W/(m K), 0.94 W/(m K), and 1.40 W/(m K) for the Clarke's models and 1.45 W/(m K), 1.40 W/(m K), 1.34 W/(m K), and 1.71 W/(m K) for the Cahill's models, respectively. Furthermore, the 3-D representation of minimum thermal conductivity of Fe3−x Cr x O4 compounds was evaluated based on the Clarke's model. Minimum thermal conductivity at (001) and (110) planes were plotted and we found that the directional dependence of thermal conductivity was determined by the anisotropy of elasticity of Fe3−x Cr x O4 compounds. This work will provide a guide for studying the property of magnetite both simulation and experiment in different Fe/Cr ratio.
ACKNOWLEDGMENTS
This work was supported by the Natural Science Foundation of Shaanxi Province of China (2013JQ6009), the 863 project in China (No. 2013AA031203), the Special Funds for Strategic Emerging Industries Core Technology Research of Guangdong Province in China (No. 2012A090100018), Postdoctoral Science Foundation funded project of China (No. 2014M552434), and Natural Science Foundation of China (51301131).







