Biaxial strain effect induced electronic structure alternation and trimeron recombination in Fe3O4

The Verwey transition in Fe3O4 is the first metal-insulator transition caused by charge ordering. However, the physical mechanism and influence factors of Verwey transition are still debated. Herewith, the strain effects on the electronic structure of low-temperature phase (LTP) Fe3O4 with P2/c and Cc symmetries are investigated by first-principles calculations. LTP Fe3O4 with each space group has a critical strain. With P2/c, Fe3O4 is sensitive to the compressive strain, but it is sensitive to tensile strain for Cc. In the critical region, the band gap of LTP Fe3O4 with both two symmetries linearly increases with strain. When strain exceeds the critical value, DOS of spin-down t2g electron at Fe(B4) with P2/c and Fe(B42) with Cc changes between dx2-y2 and dxz + dyz. The trimerons appear in Cc can be affected by strain. With a compressive strain, the correlation of trimeron along x and y axes is strengthened, but broken along the face diagonal of FeB4O4, which is opposite at the tensile strains. The results suggest that the electronic structure of Fe3O4 is tunable by strain. The narrower or wider band gap implies a lower or higher transition temperature than its bulk without strains, which also gives a glimpse of the origin of charge-orbital ordering in Fe3O4.

As an ancient magnet, Fe 3 O 4 has been used as compass with a history about 3000 years 1 . With more in-depth understanding of Fe 3 O 4 , its novel properties including half-metallicity and high Curie temperature of about 860 K have potential applications in spintronic devices 2,3 . At ambient conditions, the high-temperature phase Fe 3 O 4 has a face-center cubic lattice with a Fd m 3 symmetry. As a mixed-valence iron oxide with an inverse spinel lattice, Fe 3 4 4 . Two Fe B atoms have one spin-down t 2g electron in 3d orbits 5 . Rapid hop of the electron between two neighbor Fe B forms the conductive mechanism of Fe 3 O 4 .

O 4 is formally written as Fe
In 1939, Verwey found that the conductivity of Fe 3 O 4 drops about two orders by cooling down to 125 K 6 . The lattice symmetry transforms from cubic to monoclinic simultaneously. The first-order approximation given by Verwey is an order-disorder transition of charge distribution at Fe B 7 . The lattice structure of low-temperature phase (LTP) Fe 3 O 4 once puzzled us. Gradually, the lattice structure of LTP Fe 3 O 4 was clarified by X-ray diffraction, Raman and infrared spectroscopy in recent thirty years 1,[8][9][10][11] . Below T V , the lattice is a suppercell of × × a a a 2 2 2 c c c (a c is the cubic lattice constant) with Cc symmetry. The charge ordering has been observed by Magnetic Compton scattering 12 , resonant multiwave X-ray diffraction 13 and selected area electron diffraction 14 . With the P2/c 1,3,15 or Cc space group [16][17][18] , some theoretical calculations on the charge-orbital distribution give the results that is consistent with the experiments, which describe the ionic distribution, complex charge-orbital ordering pattern and ferroelectricity.
Recently, Senn et al. [18][19][20] proposed a new type of quasiparticle named as "trimeron" by both the experimental and theoretical results, where an anomalous shortage of the distance between Fe 2+ and Fe 3+ appears. The minority-spin t 2g electron is delocalized in a polaron that is composed of one Fe 2+ donor and two Fe 3+ acceptors [18][19][20] . Owing to the weak interactions, trimeron can be regarded as an orbital molecule, where three Fe ions locally coupled within an orbital ordered solid state. The trimeron model provides a new idea for understanding the Verwey transition. However, the case at a lattice strain may be different. High quality Fe 3 O 4 samples grown on SrTiO 3 21 and Co 2 TiO 4 22 have been investigated, where the transition temperature shows a significant difference of about 10 K. So the strain may play an important role in the transition of Fe 3 O 4 , which should be studied in details.
In perovskite oxides (ABO 3 ), B site is at the center of O octahedra 23 and covalently bonding with the nearest O atoms 24 . Some previous results show that the tilting or rotation of the O-octahedra has an influence on the band gap of perovskite 24,25 . It is well known that the change of bond angle or bond length modifies the crystal field and band structure. Borisevich et al. 25 indicate that the enlarged Fe-O-Fe angle and a higher symmetry can reduce the band gap of BiFeO 3 . By doping a ion with a larger atomic radium at B-site, Jiang et al. successfully tuned the band gap in CaFeO 3 26 .
In order to investigate the strain effect on the charge-orbital ordering and electronic structure of LTP Fe 3 O 4 , the first-principles calculations are carried out on LTP Fe 3 O 4 with P2/c and Cc. It is found that the orbital ordering and band structure of LTP Fe 3 O 4 can be tuned by the external strain. The band gap of Fe 3 O 4 with both two symmetries can be changed by a strain with a critical region, where the trimeron shows a complex relation with the external strain.

Calculational Details
The electronic structures of the LTP Fe 3 O 4 with structure (I) P2/c 1 and structure (II) Cc 19 are calculated by using the potential projector augmented wave method in Vienna Ab initio Simulation Package 27,28 . The calculations are based on the generalized gradient approximation plus on-site Column interaction (GGA + U) [29][30][31] . The energy cutoff is 400 eV. The Monkhorst-Pack grid of k-points for structure (I) is 6 × 6 × 2 and that for structure (II) is 3 × 3 × 2. The on-site Column interaction parameter U = 4.5 eV and on-site exchange interaction parameter J = 0.89 eV for all the Fe ions are used in the two structures 16 . The lattice constants and atomic positions in the two structures are used as that in refs 1 and 19, respectively. The same parameters except for k-points of 3 × 3 × 3 are used to calculate the high-temperature phase (HTP) Fe 3 O 4 with structure (III) Fd m 3 symmetry. The stress is defined by the change of lattice constants as S = (a s − a)/a × 100%, where S, a and a s represents the strain, lattice constant without and with strain, respectively. Biaxial lattice strain is applied by fixing the in-plane lattice constants (a and b) and relaxing z direction throughout the calculations. The tensile and compressive strains are defined as S > 0 or S < 0. In order to clarify the strain effects on the charge-orbital ordering, the structural optimization for structures (I) and (II) with lattice constants are carried out, where the atomic positions are fully relaxed. Then, the lattice strain is taken into considerations. The structure optimization will stop until the total energy change is less than 10 −5 eV and the Hellman-Feynman forces of optimized structure fall below 10 −2 eV/Å.

Results and Discussions
Electronic & lattice structure with P2/c symmetry. In Fig. 1, the unique equivalent site of Fe B in structure (III) breaks into Fe(B1a), Fe(B1b), Fe(B2a), Fe(B2b), Fe(B3) and Fe(B4) in structure (I) as the symmetry reduces 1,3 . The coordinate system of monoclinic structure rotates by 135° around z axis 3 . Figure 2 shows the electronic structure of structure (I) and (III). HTP Fe 3 O 4 shows a half metallic characteristic, where the spin-down states near Fermi level comes from the extra minority electron of Fe B t 2g orbits 2,5 . In Fig. 2(b), the band gap of structure (I) is opened by 0.51 eV at Fermi level, which is a bit larger than the spectroscopic 0.14 eV 10 . Figure 2(c) shows the partial DOS at different Fe B sites. The minority electron of Fe B is localized at Fe(B1a), Fe(B1b) and Fe(B4), which is consistent with previous results 1, 3,16 . These results suggest that Fe(B1a), Fe(B1b) and Fe(B4) are Fe 2+ and Fe(B2a), Fe(B2b) and Fe(B3) are Fe 3+ . In Table 1, the bond-valence sum (BVS) of Fe B is in well agreement with DOS. Herewith, the BVS expression is where R 0 is the bond-valence parameter 32 . For Fe 2+ and Fe 3+ , R 0 is 1.734 and 1.759, respectively. R i refers to the i th bond length and b is a constant of 0.37 Å 32 .
In order to investigate the strain effects, the strain of − 5%, − 2.5%, 0%, + 2.5% and + 5% are set. In Fig. 3, the orbit of spin-down t 2g electron at Fe(B4) (B4t 2g ↓ ) is almost in the xy plane at S = − 2.5%, 0%, + 2.5% and + 5%. At S = − 5%, the B4t 2g ↓ orbit shows a different style from others. In Fig. 3(a) and (b), the charge density of B4t 2g ↓ are projected onto (110) and (110) plan at S = 0% and − 5%. At S = 0%, the charge density of B4t 2g ↓ plotted on both planes shows ellipsoidal shape. At S = − 5%, the charge density of the electron plotted on (110) plane still shows ellipsoidal shape, but the charge density plotted on (110) plane shows a flower shape. This phenomena manifests  that the B4t 2g ↓ orbit lies in the (110) plane at S = − 5%. In order to figure out the critical strain, the electronic structure is also calculated at S = − 3% and − 4%. Figure 3(c) shows the B4t 2g ↓ charge density projected onto (100) plane with a strain from 0% to − 5%. The B4t 2g ↓ orbit still lies in the xy plane until the strain decreases to − 5%. Therefore, the compressive strain of − 5% is the critical value for P2/c symmetry. Meanwhile, DOS of 3d orbits of Fe(B4) also shows the same change. In Fig. 3(d), at S = − 5%, the B4t 2g ↓ orbit changes from d x 2 -y 2 to d yz + d xz by comparing the DOS at S = 0%. Actually, the B4t 2g ↓ orbit changes from d xy to d yz in HTP Fe 3 O 4 coordinate system because of the rotation of coordinate 3 .
Since the conductivity of Fe 3 O 4 is related to FeBt 2g ↓ and BO 6 distortion, it is necessary to investigate the O-octahedral distortion at different Fe B sites. In Fig. 4(a), the band gap shows a positive correlation with the increased Δ V except for S = − 5%. Herewith, Δ V is the average volume difference between Fe 2+ O 6 and Fe 3+ O 6 , E g is the energy gap near Fermi level. By linear fitting E g at different strains, we get E g = 0.946Δ V − 0.651, where E g at a compressive strain of − 5% is not included. Quantitatively, the volume of FeO 6 shows the magnitude of its crystal field. Since Fe 3+ has a stronger interaction with surrounded O 2− than Fe 2+ , the volume of Fe 3+ O 6 is smaller than that of Fe 2+ O 6 , but the electrostatic potential is higher than that of Fe 2+ O 6 . The results mean that the larger volume difference is, the more energy is needed when the electron hopes between Fe 2+ and Fe 3+ . Therefore, as the tensile strain increases, the gap becomes larger. Simultaneously, the competition between the band gap and thermal activation energy has a relation with the metal-insulator transition (MIT) temperature. Therefore, the MIT temperature of Fe 3 O 4 could be tuned by external strain. At a compressive (tensile) strain, the band gap becomes narrow (wide) and the MIT temperature of Fe 3 O 4 becomes lower (higher).
However, the above demonstration on the relation between Δ V and E g is not proper for the case at S = − 5%. Therefore, the Fe-O bond length in FeO 6 at Fe B is further analyzed. Since the charge density of Fe(B4) shows an obvious change and all the Fe(B4) atoms are equivalent with the same ambient ionic conditions, in Fig. 4  The Fe-O bond length distortion in horizontal plane also appears at Fe(B2) and Fe(B3). However, the Fe-O bond lengths at Fe(B1) do not change. Since the inversion centers and partial face centers are occupied by Fe(B1), the symmetry of Fe(B1) is higher than other Fe B sites, where the ambience of Fe(B1) is more stable than other Fe B sites. Therefore, the mode of ionic distribution does not change at Fe(B1).
In the band structure, the energy of spin-down conduction-band minimum at S = − 5% is higher by about 0.28 eV than that at S = − 4%. However, in Fig. 4(c), the valence-band maximum of Fe B t 2g ↓ is still just below Fermi level. The compressive strain of − 5% can change the structure of O-octahedra at Fe B sites. Simultaneously, the Fe-O Column interaction can raise the conduction band energy. So, in Fig. 4(a), the band gap of structure (I) at S = − 5% is larger by about 0.28 eV than that at S = − 4%.
Furthermore, the nearest six Fe-Fe distances around different Fe B sites are analyzed. Unlike Column's law, the < Fe 2+ -Fe 3+ > distance shows an anomalous shortage, which is even less than the < Fe 2+ -Fe 2+ > distance at Fe(B1) without strain. The phenomenon is consistent with trimeron model. As the tensile strain is applied, the weak bond interaction between Fe 3+ -Fe 2+ -Fe 3+ becomes tighter around Fe(B1). When the compressive strain is applied, the trimerons around Fe(B1) become weak. However, the distance between Fe(B4) and Fe(B3) becomes shorter than the Fe 2+ -Fe 2+ distance around Fe(B4). So, a more complex structure of trimeron forms, which will be demonstrated in the next section.
Electronic & lattice structure with Cc symmetry. In Fig. 5, when the symmetry reduces to Cc space group, the LTP Fe 3 O 4 lattice and the charge-orbital ordering pattern become more complex. The electronic structure of bulk without strain is firstly calculated. The band gap near Fermi level with and without structural optimization is about 1.0 and 0.7 eV, respectively. Figure 5(b) shows the total DOS of the optimized structure. The energy gap is larger than experimental result because the calculation is proceeded at 0 K. Then, the strain of − 5%, − 2.5%, + 2.5% and + 5% is applied to the Cc structure. Different from P2/c structure, the structure (II) is sensitive to tensile strain. So, we then calculated the electronic properties at S = + 3% and + 4% to figure out the critical value. Figure 5(c) shows the spin-down charge density at a height of 3c/8 and 7c/8 as a tensile strain increases from 0% to + 5%. At S < + 4%, the Fe(B42)t 2g ↓ orbit [marked with white square in Fig. 5(c)] lies in the (110) plane at  3c/8 and lies in the (110) plane at 7c/8. When the tensile strain exceeds the critical value at S = + 4% and + 5%, the Fe(B42)t 2g ↓ orbit rotates into horizontal plane. Figure 5(a) shows the atom sites of Fe(B42), which is labeled in dark blue. In Fig. 6(a), at S = + 4% and + 5% the DOS of Fe(B42) shows that the orbits of those Fe atoms change from d yz + d xz to d x 2 -y 2. The coordinate of structure (II) also rotates around z axis by 135° from cubic Fe 3 O 4 , which is consistent with structure (I). Therefore, the orbit actually changes from d xz to d xy at 3c/8, which changes from d yz to d xy at 7c/8 within cubic coordinate. In the inset of Fig. 6(a), by comparing the DOS at S = + 4% and + 5%, it is found that although the Fe(B42)t 2g ↓ orbit changes at S = + 4%, it still has the residual states projected onto d yz .
The residual states come from the out-of-plane slope of d x 2 -y 2. Then, the relationship between Δ V and E g in structure (II) is investigated. In Fig. 6, E g shows a positive relation with the increased Δ V at S < + 4%, which can be described as E g = 0.370Δ V + 0.534. However, the linear fitting parameters both slope and intercept are quite different from structure (I) due to the different structure and charge-orbital ordering pattern. Since the orbital change at Fe(B42) is obvious, the Fe-O bond length and O-octahedra distortion at Fe(B42) are taken as an example, where Fe(85) (at 3c/8) is selected as a substitute for other equivalent Fe(B42) sites.  Table 3 lists the Fe(85)-O bond lengths at different strains, revealing the reason for the orbital change at Fe(B42). The FeO 6 distorts in horizontal plane at S = + 4% and + 5%. The two shortest bonds are along (110) direction and the two longest bonds are perpendicular to (110) direction. When + 4% and + 5% strain is applied, both the shortest and longest bond are coexistent in diagonals. The Fe-O bond length along z direction suddenly decreases by about 0.1 Å when the strain increases from + 3% to + 4%. The obvious Fe-O bond length distortion in xy plane also appears at Fe(B2b1), Fe(B31), Fe(B32), Fe(B34), Fe(B41), Fe(B43) and Fe(B44). Due to the tensile strain, the expansion of the equatorial plane of O-octahedra can release the electrostatic energy between the surrounded O 2− and outside electron of Fe 2+ . Correspondingly, the Column interaction along z direction can also be weakened by the transformation of Fe(B42)t 2g ↓ orbit, so the Fe-O bond length along z direction becomes shorter. Since the equatorial face can further expansion at S = + 5%, more electrostatic energy can be released, where the Fe(B42)t 2g ↓ orbit becomes more parallel to the xy plane. In the inset of Fig. 6(a), at S = + 5%, the residual d yz states are less than that at S = + 4%. In Fig. 7(b), the conduction-band minimum at S = + 4% is lower than that at S = + 3%, where the valence-band maximum is still just below Fermi level. As a result, the band gap becomes smaller as the strain increases. Furthermore, the model of trimeron presented by Senn et al. [18][19][20] is also observed in our calculations. It is found that the distribution of trimeron can be affected by external strain. When the strain increases from 0% to + 5%, the Fe-Fe distance along x and y axis changes faster than that along the face diagonal direction of Fe B4 O 4 . The Fe-Fe distance along diagonal direction even reduces with the increased strain at some Fe B sites. The process of Fe B4 O 4 distortion is compared with an ideal model of equivalent volume deformation in tetragonal system. . At 0 < a < 3.367, < 0 dl da . In our model, the case of 0 < a < 3 is considered. As a result, the length of face diagonal reduces with the increased lattice constant. So, the   trimerons along x and y direction break down by the tensile strain, but the correlation of trimerons along the face diagonal are strengthened. When a compressive strain is applied, the Fe-Fe distance along x and y direction becomes short and the Fe-Fe distance along face diagonal elongates. The trimerons along x and y direction are strengthened, but the trimerons along the face diagonal directions break down due to the compressive strain.

Conclusions
We have investigated the biaxial strain effects on the electronic structure of LTP Fe 3 O 4 with P2/c and Cc space group by GGA + U method. When the strain on the two structures are below their critical region, the distortion of O-octahedra can change the electrical potential difference between the nearest ferric and ferrous ions. As a result, the band gap shows a positive linear correlation with the strain. The narrower or wider band gap implies a lower or higher transition temperature. When the strain is above the critical value namely S < − 4% in structure (I), the orbit of Fe(B4)t 2g ↓ changes from d xy to d yz in HTP Fe 3 O 4 coordinate and the energy of conduction-band minimum raises. In structure (II), at S ≥ + 4%, the orbit of Fe(B42)t 2g ↓ changes from d yz to d xy in HTP Fe 3 O 4 coordinate and the energy of conduction-band minimum reduces. The trimeron appears in both the structure (I) and (II). The distribution of trimeron can also be affected by strain. The trimerons along x and y axes get broken (strengthen) at a tensile (compressive) strain. However, the trimerons along face diagonal are broken (strengthened) at a compressive (tensile) strain. These results can be ascribed to the change of Fe-Fe distance when different strains are applied, which can be estimated by geometric calculations.