Magnetic ordering induced giant optical property change in tetragonal BiFeO3

Magnetic ordering could have significant influence on band structures, spin-dependent transport, and other important properties of materials. Its measurement, especially for the case of antiferromagnetic (AFM) ordering, however, is generally difficult to be achieved. Here we demonstrate the feasibility of magnetic ordering detection using a noncontact and nondestructive optical method. Taking the tetragonal BiFeO3 (BFO) as an example and combining density functional theory calculations with tight-binding models, we find that when BFO changes from C1-type to G-type AFM phase, the top of valance band shifts from the Z point to Γ point, which makes the original direct band gap become indirect. This can be explained by Slater-Koster parameters using the Harrison approach. The impact of magnetic ordering on band dispersion dramatically changes the optical properties. For the linear ones, the energy shift of the optical band gap could be as large as 0.4 eV. As for the nonlinear ones, the change is even larger. The second-harmonic generation coefficient d33 of G-AFM becomes more than 13 times smaller than that of C1-AFM case. Finally, we propose a practical way to distinguish the two AFM phases of BFO using the optical method, which is of great importance in next-generation information storage technologies.

Magnetic materials have received intensive attention [1][2][3][4][5][6][7][8][9][10][11][12] in the last decades, as they generally contain rich mutual interactions of spin, charge, photon and lattice degree of freedom. They also have important applications in advanced technologies such as magnetic read heads and magnetic memory cells, and many other new tools for magnetic imaging [13][14][15][16] . Despite extensive studies in the field of magnetic materials, the tools to obtain detailed information on the arrangement of magnetic moments in antiferromagnets are still limited, largely due to the fact that antiferromagnets have no or little macroscopic magnetization and therefore are insensitive to the external magnetic field 17,18 . Considering the magnetic ordering often has important effect on the properties of the material [19][20][21][22][23][24][25][26] , a practical, inexpensive way to its measurement is in urgent demand, which is also important to both fundamental and technological developments in the field of magnetoelectronics and spintronics.
As an invaluable experimental probe, the optical spectroscopy is of great scientific and practical interest due to its advantage in providing complementary information on crystallographic, electronic and magnetic properties and studying the coexistence and interactions of magnetic and electric order 27,28 . With the rapid development of spintronics and its cousin valleytronics, the optical spectroscopy is regarded as a promising candidate to probe and manipulate the spin and valley degrees of freedom [29][30][31][32][33] . Compared with neutron scattering method, which is regarded as the most powerful and versatile experimental technique in obtaining truly three-dimensional pictures of magnetic structures in various materials 34 , optical spectroscopy techniques possess the advantage of relative small samples, simplicity of source preparation, rapidity of measurement in a non-contact way. Since the arrangement of spins has great influence on the interplay between matter and electromagnetic radiation, which will directly be reflected in optical properties of materials, is it possible to determine the magnetic ordering using optical methods?
In this work, we focus our study on the perovskite oxide bismuth ferrite (BiFeO 3 , BFO), which is the very rare single-phase multiferroic material simultaneously with antiferromagnetism and robust ferroelectricity well above room temperature 35,36 , and has attracted extensive research [37][38][39][40][41][42][43][44][45][46][47][48][49] . Using first-principles calculations, we investigate the complex dielectric functions and second-harmonic generation (SHG) coefficients of expitaxially stabilized BFO in tetragonal phase with giant c/a ratio and find sizable differences for both linear and nonlinear optical properties between C1-and G-antiferromagnetic (AFM) ordering. The predicted modification of optical properties induced by magnetic ordering change reveals the possibility of using a noncontact and nondestructive optical method to distinguish the C1-and G-type AFM of BFO, which is currently hard to be handled by a traditional neutron diffraction measurement because of the strong diamagnetic response of substrate materials 50 .

Results and Discussion
A tetragonal BFO with giant c/a ratio, i.e. the so-called super-tetragonal structure 51 , is more energetically favorable under compressive strains with an "out of plane" ferroelectric polarization 41 . There exists a transition between C1-and G-type AFM ordering for tetragonal BFO under the influence of in-plane strain 19,52 . As a result, here we establish the [001] polarized BFO in tetragonal perovskite structure (C 4v , P4 mm) with the lattice constant of 3.905 Å, around which is the transformation point for the two AFM phases. The optimized c/a ratio is 1.197.
To explore the effect of magnetic ordering on the electronic structure, the spin-dependent charge density is needed to be inspected. As an example, in Fig. 1, we display the majority-spin ones. Noted that on account of AFM ordering for both C1 and G states, the minority-spin charge density plots can be directly obtained by shifting the corresponding majority-spin ones by the in-plane lattice constant a along the x-or equivalently y-axis. Since all the eight operations in point group C 4v are z-irrelevant, the same intraplane couplings of C1 and G cases signify the same magnetic point group, i.e. C 2v . In xy plane, the magnetic coupling of neighboring magnetic ions in both states are AFM, which leads to essentially the same charge density distribution. On the contrary, when we consider the z direction, there exists significant difference in the spin-dependent charge density of Fe atoms. In C1-AFM, whose nearest interlayer coupling is ferromagnetic (FM), as shown in Fig. 1(c), the density changes continuously along the z direction in the spin-up channel, and so does the spin-down one. For G-AFM, however, the nearest AFM interlayer coupling leads to the discontinuous charge density distribution, which inevitably affect the spin-transport properties of the BFO film.
The spin arrangement along the z-direction not only have influence on transport behaviors, but also will affect the band structure. The energy spectra E n (k) near the Fermi level are displayed in Fig. 2. We emphasize that to make a straightforward comparison between the band structures of two AFM phases, the density functional theory (DFT) calculated bands along high-symmetry directions are identical in the k-space, although their first Brillouin zones (FBZ) are indeed different. In both C1-and G-type AFM orderings, the ingredients of electronic states are quite similar. The presented valance band (VB) close to the Fermi level is predominantly comprised of Fe-e g and O-2p states. The band mainly from the Bi-p z orbital is highly dispersive along Γ -Z and Z-A directions. The flat band lying around 1.7 eV above the Fermi level is primarily of the empty Fe-3d xy character. Nevertheless, the dispersions of energy bands of the two magnetic cases are significantly different. For C1-type AFM, the top edge of VB, as well as the bottom of conduction band (CB) are located at the Z (0, 0, 0.5) point. The band gap (E g ) is direct, about 1.12 eV. For G-type AFM, although the lowest unoccupied molecular orbital (LUMO) is still at the Z point, the highest occupied molecular orbital (HOMO) shifts to the center of the FBZ. Its E g is indirect with the value of 1.03 eV. The direct gap E g is much larger, i.e. 1.51 eV at the Z point. When we carefully analyze the band features, especially along the Γ -Z symmetry line, it is intriguing to point out that due to the occupancy of the Fe-3d states, the extreme points of bands would be shifted during the transformation of magnetic ordering. As an example, for C1-AFM, the maximum and minimum of the d z 2 and d xy bands are located at Z and Γ , respectively. For G-AFM, they are just reversed. However, it would not happen for the band dominated by the p electrons of Bi atoms. In order to accurately capture the physics about the impact of Fe-3d orbitals on the energy band dispersions between two AFM states, particularly the shift of the HOMO, tight-binding (TB) model calculations are then carried out.
As our purpose is to show the band difference between different magnetic configurations, we adopt an effective, iron-only one-band TB model. This is based on the fact that except for the mixing between d z 2 and − d x y 2 2 orbitals along Γ -A direction in BFO under G-AFM ordering, the d xy , d z 2 and − d x y 2 2 states are basically independent, as shown in Fig. 2. Since the electrons hopping between Fe ions generally needs the O-2p orbitals as medium, the iron-oxygen hybridization has already been included in the effective hopping term of our iron-only model. Particularly, the minimal TB model involving Fed z 2 orbital is used to reveal the behavior of the top edge of VB.
To be specific, the interatomic matrix elements of the non-magnetic Hamiltonian can be written in terms of on-site energies ε aα and effective hopping integrals t bβ,aα as follows: Here k is the Bloch wave vector. a (b) and α (β) are the atom index and the atomic orbital quantum number, respectively. In order to take the magnetic coupling into account, Hilbert space should be doubled by including the spin degree of freedom (σ), i.e. |aαk> → |aασk> . Then the effective hopping integral between any two ions depends on cos(θ/2) 53 , where θ is the angle between their spins. Consequently, i.e. only when the coupling is FM, there exists hopping between the two magnetic ions and the interaction of them should then be invoked in the Hamiltonian. Considering that for tetragonal BFO under C1-and G-type AFM ordering, the magnetic coupling of Fe ions in the xy plane are exactly the same, we write the TB Hamiltonian referring to the Fe-3d orbitals as: and R i is the lattice constant along the x, y or z direction. ε d refers to the single particle energy. We consider the intralayer and interlayer hopping between first, second and third FM neighbors for C1-and G-AFM, accordingly. Test calculations show that higher order hopping terms are negligible and therefore are not counted. The interaction pathways corresponding to the effective hopping parameters are shown in Fig. 3.
A straightforward analysis show that for both magnetic configuration the band reaches its extreme point at Γ and Z. The energy differences between them are given by: Obviously, locations of the maximum and minimum for the band are determined by the interlayer hopping parameters.
As listed in Table 1, effective hopping integrals of the one-band TB model involving Fed z 2 orbitals are obtained by fitting to the first-principles results. For C1-type AFM case, the negative parameter t 1c ′ causes that the subband has higher energy in the Z point than in the Γ point. For G-AFM, however, the interlayer nearest hopping parameter t 1c ″ is positive and then leads to the opposite results. The two-center Slater-Koster parameters 54 using the Harrison approach 55 are applied to explain the sign difference between them. For d z 2 orbitals, the hopping term is a function of the direction cosines, l, m, n of the vector x, y and z, as shown below. The three dd matrix elements are all obtained from   where r d is a length that is characteristic of each element and d is the distance between atoms. ћ and m are reduced Planck constant and electron mass, respectively. Then the terms t 1c ′ and t 1c ″ can be written as: For the term t 1c ′ , since the interaction path is just along the z direction, only σ-orbital combination exists between the Fe atomic orbitals. Same signs of the wave functions makes the coefficient η ddσ negative, and then gives rise to the negative t 1c ′ . While, for the term t 1c ″ , there are angular momentum around all the three vectors, which triggers the coexistence of σ, π and δ orbital combinations. The coefficient η ddδ is negligible compared with η ddσ and η ddπ . The π-orbital combination has factor 3(c/a) 2 , which is almost 5 times larger than the factor ((c/a) 2 − 1/2) 2 for σ-orbital combination. Therefore, the σ-orbital combination is suppressed by the π-orbital combination. The positive coefficient η ddπ makes the term t 1c ″ positive. We also examine other effective hopping parameters in Table 1. The signs of them are all consistent with that of the two-center Slater-Koster parameters using the Harrison approach.
We emphasize that the effective hopping parameters t 1c ′ and t 1c ″ are the most important fitting parameters to explain the behavior of the d z 2 band. Sign difference between them directly controls the location of the HOMO. When the magnetic ordering of tetragonal BFO changes from C1-to G-type, the interlayer nearest hopping term varies from the negative t 1c ′ to the positive t 1c ″ , directly triggering the top of VB shifts from the Z point to the center of the FBZ.
In addition to the d z 2 orbital we discussed above, the minimal TB model involving − d x y 2 2 and d xy orbitals of Fe-3d are also taken into account. Figure 4 compares the model TB bands (open circles) with the corresponding DFT results (solid lines). Due to the hybridization between d z 2 and − d x y 2 2 states, which is not considered in our simple one-band TB models, there exists mismatch in bands dispersion along the Z-A symmetry line for G-type BFO. Even so, the agreement between the TB models and the first principles calculations is excellent, which warrants the validity of the TB models we adopted here.
According to previous discussions, the influence of magnetic ordering on band structures primarily reflects on the shift of the HOMO, the reason of which has been explained by the simple yet convincing TB models. As a result, in the process of phase transformation from C1-to G-AFM in tetragonal BFO, the direct band gap becomes indirect. We notice that the nature of the band gap of BFO both in tetragonal and rhombohedral phase 38,46,47 remains debatable. In general, temperature, disorder or strain effect could dramatically affect whether the band gap is direct or not. Here, by combining the first-principles calculations with TB models, we propose a new mechanism, i.e. magnetic ordering could also have significant influence on the band gap and even the whole band structures. Since the magnetic orders in experimental work are generally moot, we could not assert but prefer to indicate that the AFM ordering transition stands a good chance to cause the controversy.
It is interesting to point out that, when the magnetic ordering of tetragonal BFO varies from C1-to G-type, in contrast to the decrease of band gap (~0.09 eV), the enhancement of direct band gap reaches up to 0.40 eV. It, apparently, will greatly affect the optical properties of BFO. For the linear optical properties, the imaginary parts of complex dielectric functions ε 2 are calculated and shown in Fig. 5. Note that as a second-order process with a relatively small transition probability, the indirect transition process would happen only if phonons participate in and is therefore ignored here. The interband transitions begin at 1.12 eV and 1.52 eV for C1-AFM and G-AFM, respectively, in good agreement with the direct E g gained by band structures. For C1-type AFM case, owing to the selection rules, the transition from HOMO to LOMO in the Z point is allowed when the incident radiation is polarized in z-direction. The optical band gap (E g opt ) exactly equals to the direct E g . However, for G-type AFM case, for any polarization direction of incident light, the electric-dipole transition related to the direct E g in the Z point is forbidden. Nevertheless, it is allowed when the k-point moves to the Z-A symmetry line. This suggests that E g opt in this case should be slightly larger than the already increased direct E g . Consequently, as shown clearly in Fig. 5, when the magnetic ordering varies from C1-to G-type AFM, E g opt of tetragonal BFO exhibits a blue shift of approximately 0.4 ~ 0.5 eV, which is large enough to be easily distinguished in experimental optical spectra.
Besides the complex dielectric functions, the nonlinear optical properties, especially in the static limit, are more sensitive to the change of the direct E g . The SHG susceptibilities in static limit are listed in Table 2. Due to the symmetry of tetragonal BFO with the point group of C 4v , there are five independent SHG components, i.e. d 15 = d 24 , d 31 = d 32 and d 33 . The so-called Kleinman symmetry 56 demands d 15 = d 31 in the static limit. Though the d 15 component between the two kinds of AFM ordering is about the same, the absolute value of d 33 for C1-type AFM phase is more than 13 times larger than that of G-type AFM phase. In addition, what is noteworthy is that there exists considerable difference between the component d 15 and d 33 with an order of magnitude in G-type BFO. While the similar phenomenon cannot be seen in the case of C1-type AFM ordering. Such evident optical anisotropy in SHG coefficients in static limit reveals its utility to measure the magnetic ordering in experiments by using different polarization combination to determine the SHG procession efficiency in terahertz (THz) of the tetragonal BFO.
Finally we propose a practical device to distinguish the C1-and G-type AFM of BFO using optical methods. As shown in Fig. 6, the tetragonal BFO under the [001] polarization direction with the lattice constant 3.905 Å, around which the C1-G magnetic transition will occur, is placed onto a PMN-PT ([Pb(Mg 1/3 Nb 2/3 )O 3 ] 0.72 [PbTiO 3 ] 0.28 ) actuator via gluelike PMMA (polymethyl methacrylate) layer. The biaxial strain is provided by the PMN-PT substrate sandwiched between transparent In 2 O 3 thin films acting as electrodes. A bias voltage V applied to the PMN-PT results in an out-of-plane electric field E, which leads to an in-plane strain ε. The PMN-PT is electrically poled so that V > 0 (< 0) corresponds to in-plane compressive (tensile) strain. Several experimental works 57,58 demonstrated the reversibility of the strain tuning technique. In order not to affect the transmission of light, the materials serving as actuator, glue and electrode are carefully selected to be highly near-infrared transparent. A laser beam with the energy hv = 1.3 eV is incident to the tetragonal BFO. Note that as DFT calculation generally will underestimate the energy gap, the realistic photon energy could be larger. For V > 0, the tetragonal BFO experiences an in-plane compression and is in the C1-AFM state. Under this circumstance, electrons will absorb the photons and jump to the excited state, i.e. the film is opaque. By applying a negative voltage V < 0 to the substrate, however, the BFO film transforms to G-type AFM. Since now the energy gap is larger than the photon energy, the film then becomes transparent, and the incident light could be detected by the IR detector beneath.  Table 2. SHG coefficients d ij (pm/V) in the static limit for the tetragonal BFO with C1-and G-AFM ordering.
In the above device, the magnetic orderings can be controlled by the sign of the bias voltage V. And the magnetic states can be easily "read out" utilizing the IR detector to monitor the transmitted light. Obviously, it might also be used as the electrically writing and optically reading memory devices.

Conclusion
Using first-principles calculations, we found that the different spin arrangements in tetragonal BFO have significant influence on the dispersion of energy bands, especially the location of the HOMO, which has been explained by the minimal one-band TB models. When the magnetic ordering of tetragonal BFO varies from C1-to G-type, sign difference between the interlayer nearest hopping terms directly triggers the shift of HOMO from the Z point to Γ point. Consequently, the energy gap changes from direct to indirect one and leads to drastic modification to the optical properties of tetragonal BFO. For the linear ones, the enhancement of E g opt can be as high as 0.4 eV. When it comes to the nonlinear ones, the change of the SHG coefficient d 33 becomes more than 13 times. Obviously, the difference of optical properties are considerable enough to distinguish the C1-type and G-type AFM of tetragonal BFO experimentally. We therefore hope that our theoretical predictions will stimulate experimental studies of magnetic ordering measurements using a noncontact and nondestructive optical method.

Methods
The calculations are performed within DFT using the accurate full-potential projector augmented wave (PAW) method, as implemented in the Vienna ab initio Simulation Package (VASP) 59 . The exchange-correlation potential is treated in Perdew-Burke-Ernzerhof (PBE) form 60 of the generalized gradient approximation (GGA) with a kinetic-energy cutoff of 500 eV. A 7 × 7 × 7 and 11 × 11 × 11 Monkhorst-Pack k-point mesh centered at Γ are respectively adopted in the self-consistent and optical calculations. A primitive cell, containing two BFO formula units (10 atoms) is performed to describe the magnetic structures of BFO. Due to the different shapes of the primitive lattice between C1-and G-type AFM cases, all the coordinates mentioned below are described in real or reciprocal space of C1-AFM. The convergence criterion for the electronic energy is 10 −6 eV and the structures are relaxed until the Hellmann-Feynman forces on each atoms become less than 1 meV/Å. We utilize the LSDA + U scheme 61 , which is adequate for searching magnetically ordered insulating ground states 62,63 by taking the effect of Coulomb repulsion into account and the effective Hubbard constant U eff = U -J = 2.0 eV 49 is chosen to give a better description of the electron-electron interaction between the partially filled and strongly correlated localized Fe-3d electrons. We also check that our results are robust within the other widely used U eff 's as 4.0 eV and 6.0 eV.
The optical calculations are performed using our own code OPTICPACK, which has been successfully applied to study the optical properties of borate crystals 64,65 and other multiferroic materials 66 . For the linear case, the imaginary part of the complex dielectric function ε 2 is calculated using the following relations: where E is the photon energy, Ω is the cell volume, p cv is the electron momentum matrix element between the VB states (v) and the CB states (c). The integral over the k space has been replaced by a summation over special k points with corresponding weighting factor  W k . The second summation includes v and c states, based on the reasonable assumption that the VB is fully occupied, while the CB is empty.
For the nonlinear case, we only consider the SHG susceptibility χ (2) (2ω; ω, ω) in the static limit, which can be written as:  in which P(ijk) denotes full permutation over Cartesian components i, j, k and explicitly shows Kleinman symmetry 56 of the SHG coefficients. The first and second summations represent virtual-hole and virtual-electron process respectively.