A Dixit

D Stoeffler

M Alouani

Effect of site disorder on the Electronic, Magnetic and ferroelectric Properties of Gallium Ferrite

The electronic, and magnetic properties of the magnetoelectric solid GaFeO 3 are investigated within the generalized gradient approximation including the Hubbard interaction (U) on the localized d orbitals of iron. It was found that using an on-site U =8 eV describes consistently the experimental results. The origin of ferrimagnetism was attributed to the cationic site disorders.

The density of states at the iron sites, in a octahedral geometry, show that the occupied e g states are below the t 2g states in contradiction with the crystal-field splitting obtained by a point charge model. However for the unoccupied states the ab initio data agree qualitatively with the model, showing the complexity of the electron-electron interaction in GaFeO 3 . The computed electric polarization of the system as a function of the temperature in the linear regime shows a monotonic decreasing trend. To determine the nature of the magnetoelectric coupling, the computed electric polarization as a function of the rotation of the magnetization axis indicate that the magnetoelectric effect observed experimentally could not have been due to a direct coupling between the electric and magnetic order parameters. Finally the calculated x-ray absorption and x-ray magnetic circular dichroism spectra for the disordered system is shown to be in good agreement with experiment.

I. INTRODUCTION

Magnetoelectric materials are multiferroic materials, which have coupled magnetic and electric ferroic order parameters. This type of coupling in materials was speculated by Curie [1,3] as early as in 1894. However, due to the difficulty of combining magnetic and electric ferroic orders in materials, this field was not pursued further. The first theoretical prediction of the magnetoelectric coupling was done by Dzyaloshinskii [4] in 1959 for chromium oxide, and was soon observed by Astrov [5] in 1960.

The cross-play of the ferroic properties has motivated further research into potential technological applications, especially those where the electric properties are controlled with magnetic fields or vice versa [6]. Few materials have been reported to show a sizeable magnetoelectric effect, among which gallium ferrite (Ga 2-x Fe x O 3 or GFO) appears to be of considerable significance due to the coexistence and coupling of a magnetization and electric polarization at room temperature [7].

The first gallium ferrite crystals, Ga 2-x Fe x O 3 with 0.7 ≤ x ≤ 1.4, were synthesized by Remeika in 1959, and were described as a ferromagnetic-piezoelectric material [8]. The structural characterization of GFO was determined by Wood, who found the crystallographic space group to be P c2 1 n [9]. This was confirmed by Abrahams et al. in 1965, with the lattice constants a = 8.7512 ± 0.00008 Å, b = 9.3993 ± 0.00003 Å and c = 5.0806 ± 0.0002 Å [10]. There are four different cationic sites occupied by Fe and Ga cations: three irregular (distorted) octahedral sites (Fe1, Fe2 and Ga2) and a regular tetrahedral site (Ga1) oriented along the b-axis. The oxygen anions are positioned in six different sites in a double hexagonal compact arrangement (Fig. 1).

Other studies [11][12][13] were carried out to probe the magnetic structure and magnetoelectric characteristics of GFO. However, it was in 1965 that Frankel et al. using Mössbauer spectroscopy, showed a ferrimagnetic order with magnetic moments aligned almost along the c-axis, instead of a canted antiferromagnetic order [14]. Both cationic sites Ga1 and Fe1 are antiferromagnetically coupled to the Ga2 and Fe2 sites, which should result in a net antiferromagnetic configuration for Fe composition x = 1. The presence of ferrimagnetism without the presence of a canted antiferromagnetic configuration hints at a possible existence of site disorder. In 2004, Arima et al. found that the GFO preparation method influenced the magnetic transition temperature, which were likely caused by the change in the Ga/Fe occupations at the four cationic sites [15]. The origin of ferrimagnetism was interpreted as follows -the magnetic Fe cations at the Ga1 and Fe1 sites adopt a magnetic orientation antiparallel to the magnetic cations at the Ga2 and Fe2 sites through a superexchange mechanism mediated through the oxygen anions as described in Ref. [16]. Since, the amount of Fe at the Fe2 and Ga2 sites is larger than that at the Fe1 and Ga1 sites, there is a net non-zero magnetic moment along the c-axis (parallel to the magnetic moment on Fe2). They also demonstrated experimentally that the magnetic transition temperature could be tuned with the Fe/Ga ratio and it increased with Fe content to room temperature for x ≥ 1.1. In addition, a large linear magnetoelectric effect was measured for GFO single crystals, which was one order of magnitude larger than the value reported for Cr 2 O 3 .

In 2006, Kim et al. reported a large orbital moment in GFO [17]. This was unusual because Fe in GaFeO 3 has a formal valence of +3 (half-filled d 5 configuration), for which the orbital moment is expected to be zero. The mechanism behind the orbital moment could help provide a better understanding of the magnetoelectric coupling in this material.

Despite several experimental studies [11][12][13][14][15]17], very few theoretical reports investigating the properties of GFO are present in literature [16,[18][19][20][21]. It is therefore important to ask what is the most appropriate method for the determination of the electronic structure of a large band-gap magnetic insulator? There are of course many accurate methods, such as the dynamical meanfield theory(DMFT) [22], or DMFT starting form a GW calculation, the socalled DMFT+GW [23], but those methods are computationally prohibitive for materials with many atoms per unit cell. One of the most interesting method, which provides quantitative results, without the high CPU cost, is the so called LDA+U method [24]. Indeed, it has been shown that this method is capable of providing accurately the electronic properties of materials with localized orbitals, such as transition-metal oxides, without high computational costs [24][25][26]. The only issue with this method, is how to choose the Hubbard U parameter for a material? Anisomov and coworkers [24][25][26] have shown that the U parameter can be determined for each material by constrained DFT calculation. In general, for oxides the values of U are shown to be large, but nevertheless can reproduce photoemission, BIS [25], valence bands of oxygen Kα x-ray emission, and X-ray photoemission [26]. Other calculations have shown that for CoO not only the band gap and the spin and orbital moments are reproduced , but also the magnetic-anisotropy energy and the orientation of the magnetic moment of CoO under strain caused by silver or MnO substrates [27]. However, Rödl and coworkers [28] claim that may be smaller values of U should be used for oxides. They showed that for a smaller value of U , the LDA+U is in agreement with a GW calculation when starting from hybrid HSE functional (GW+HSE), provided that the LDA+U band gap is increased by hand (the so-called scissor operator) to achieve agreement with GW+HSE results. However, the latter results are not in better agreement with experiment compared with the LDA+U results for higher values of U in the literature. One can indeed check that the results of photoemission of NiO of Refs. 25 and 26 are in much better agreement with experiment than those of Ref. 28. Since we don't want to adjust the band gap of GFO by hand, we have chosen in this paper to follow the prescription of Anisimov et coworkers and use a higher value of U to describe the band gap and the magnetic properties of GFO. In addition, our previous LDA+U results [18] with higher value of U described correctly the band gaps and the magnetic properties of GFO with excess iron up to 40%. Nevertheless, the results of Rödl and coworkers are interesting, and might suggest to combine the GW method and LDA+U to describe better the short and long range exchange and correlation potential as already performed in the GW+DMFT method [23].

Early density-functional theory (DFT) calculations carried out on the ideal structure revealed a stable antiferromagnetic state with zero net spin and orbital moments [29]. Han et al. also showed that the energy difference between an ideal GFO structure and a structure with an Fe interchanged with the Ga2 site can be as small as 1 meV per unit cell, thereby implying that this kind of site disorder is highly probable and in accord with the presence of Fe at the Ga2 sites reported in experiments [29]. With the help of first-principle calculations, Roy et al. showed that the site disorders are not favored in the ground state and that available thermal energy at room temperature (kT ∼ 25 meV) is of the order of the energy difference for the Fe2-Ga2 site disorder, hinting towards the role of thermally induced defects [30]. The same group calculated the electronic structure and the Born effective charges, which showed a largely ionic character of the Ga/Fe-O bonds and a lack of significant anomaly in the Born effective charges [31]. In 2012, Stoeffler calculated the electric polarization of the system to be -25µC/cm 2 , an order of magnitude larger than the value estimated by Arima et al. [15,32] along the b-axis. In the same year, Hatnean et al. reported a weak dependence of the Ga/Fe disorder and the magnetic transition temperature on the growth conditions, in contrast to a strong dependence reported by Arima et al. [15,33]. They also reported that the disorder affected magnetic excitations, evidenced through the damped spin waves.

To elucidate the above magnetoelectric properties of bulk GFO, we performed ab initio studies, under different approximations, including spin-orbit coupling, to probe the electronic, magnetic and magnetoelectric properties for the iron concentration, x = 1. Our analyze of the electronic properties for the structure with cationic site disorders provided the origin of ferrimagnetism and a good comparison with the experimental results by Hatnean et al. [33]. Our investigation of the magnetoelectric effects in the ideal structure by changing the direction of the magnetic field shows that the coupling between magnetization and the electric polarization is weak contrary to the experimental results. Finally the crystal-field splitting is understood and analyzed in terms of our ab initio calculations and a point charge model. Our calculations for the x-ray absorption spectra (XAS) and x-ray magnetic circular dichroism (XMCD) spectra for a disordered system for comparison with the spectra obtained by Kim et al. [17] showing that the disordered structure used is correct.

Our paper is organized as follows, in the second section we briefly describe the details of our DFT calculations used to obtain the electronic and magnetic properties of GFO. In the third section we give first the electronic and magnetic properties of ideal GFO and then the same properties with site disorder to determine the origin of the ferrimagnetism and at the end of this section we address the crystal-field splitting and compare our ab initio results to a point charge model. In the fourth section we analyze the magnetoelectric properties of GFO and their temperature dependence in the ideal structure by changing the direction of the magnetic field. In the fifth section we present our calculations for the XAS and XMCD spectra for a disordered system and compare them to the experimental spectra of Kim et al. [17].

II. METHOD OF CALCULATIONS

Our DFT calculations were carried out using the VASP package [34,[START_REF]The vasp manual[END_REF]. We used the projector augmented wave basis set [START_REF] Blöchl | [END_REF]37] and the exchange-correlation functional was described using the local density approximation (LDA) as parameterized by Perdew and Zunger [38], and the generalized gradient approximation (GGA) as parameterized by Perdew, Becke and Ernzerhof [39,40]. We used k-point mesh of 7 × 7 × 9 in the Brillouin zone, which is required to converge the ground state energy and the magnetocrystalline anisotropy energy to within 10 µeV [18]. For the plane wave cut-off, a value of 550 eV was used and the convergence criterion for the electronic self-consistent loop for the total energy was set to 10 µeV. Both LDA and GGA underestimate the energy band gap and magnetic moments of GFO as compared to the experimental values. This failure of LDA and GGA is known in correlated systems like transition metal oxides [41][42][43]. To accurately account for the strong on-site Coulomb interaction among the localized Fe 3d electrons, we used the rotationally invariant approach introduced by Lichtenstein et al., represented by the Hubbard-like term U and the exchange term J [44]. This led to an improvement of the ground state properties of GFO. Based on an earlier theoretical study, the value of J was set to 0.9 eV [18]. In this work, two specific values of U were used, U = 4 and 8 eV, chosen so as to probe the ground state properties as well as the hybridization of the localized orbitals with the remaining delocalized states. The spin-orbit coupling (SOC) was included in some of our calculations as implemented in VASP to improve the magnetic properties of the system [45] and to calculate the x-ray magnetic circular dichroism.

To study the origin of ferrimagnetism in GFO, the ionic occupancies provided by Hatnean et al. were used [33]. In addition, a simple predictive model was developed that allowed us to compute the net spin and orbital magnetic moments in a disordered structure with the help of the cationic occupancies. The predictions of this model were then compared with the values obtained by the ab initio calculations and experiment.

For the crystal-field analysis, the global frame , m-site-projected basis set were rotated to the local octahedral or tetrahedral frame of reference to better understand the symmetries of the e g and t 2g sub-orbitals. The rotation was first defined directly using the Euler angles and the real spherical harmonics (also known as the cubic harmonics) and then implemented in the VASP code. To study the effect of crystal-field and hybridization on the e g -t 2g splitting on the octahedral Fe sites, both values of U were used. The e g -t 2g splitting was modeled using a point-charge model, where the potential of the oxygen ligands was expanded in terms of spherical harmonics in the octahedral center. Such a model allowed us to determine the effect of the ligand on the splitting of iron 3d states. The diagonalization of this matrix, using the cubic harmonics basis set, gave us the e g -t 2g splitting which we discussed and compared to our ab initio calculations.

The electric polarization of GFO was computed as implemented in the VASP code [46][47][48][49][START_REF] Resta | Berry Phase in Electronic Wavefunctions[END_REF][51].

The preliminary calculations were done with the aim to study the temperature dependent electric polarization on a path between the experimental positions at 4 K and 230 K provided by Arima et al. [15], and to recover the polarization calculations done previously [32]. This initial calculations were extended to probe the magnetoelectric effect, where the magnetization axis was rotated for different angles along the c-b plane.

III. ELECTRONIC AND MAGNETIC PROPERTIES OF GFO A. Ideal bulk GFO

In this section, we present the electronic and magnetic properties of ideal GFO. The lattice parameters and the atomic positions were based on the values reported by Arima et al. [15], which were obtained with neutron diffraction studies at 4 K. The values obtained at 230 K were used in the later magnetoelectric studies. The lattice parameters are a = 8.719 Å, b = 9.368 Å and c = 5.067 Å. The atoms were not relaxed since the parameters were obtained from experiment. Previous calculations have shown that atomic relaxations do not change significantly the experimental positions [18]. Two values of U were used for the Fe 3d orbitals: 4 eV and 8 eV, for both LDA and GGA, to improve the ground state properties like the energy band gap and magnetic moments. The calculations including the SOC does not significantly affect the physical properties of interest, but allows us to compute the orbital magnetic moments and XMCD and are of importance for the magnetoelectric effects.

To understand the electronic structure, we plot the total electronic density of states (DOS) in Fig. 2. When the value of U under both LDA and GGA is increased, the Fe 3d orbitals become more localized and the upper and lower Hubbard bands are split further. This is in accordance with the Hubbard model, where higher values of U tend to localize the d orbitals. Increasing the U value beyond 8 eV will eventually produce the electronic structure of an isolated Fe atom. Also, it should be noted that under both LDA and GGA, though the value of U = 8 eV gives good values of the energy band gap and the spin magnetic moments, the hybridization of the Fe 3d orbitals with the 2p orbitals of the neighboring O atoms is drastically reduced and may not be physical. Since no photoemission data are available to compare with our DOS, we used the value of U that best describes the observed properties like the energy band gap and the magnetic moments, i.e. U = 8 eV. The values for these properties are tabulated in Table I.

It can be seen that value of U = 8 eV for LDA+U and GGA+U yields an energy band gap in good agreement with experiment and an improvement in the prediction of the magnetic moment for the Fe2 site. However, the magnetic moment at the Fe1 site does not agree with experiment. As is shown in the subsequent section, this discrepancy arises due to the partial iron occupancy at the Fe1 site, which was estimated experimentally to be 84%. Assuming a spin moment at 84% of -4.44 µ B gives us a value of -3.73 µ B , which is closer to experiment. In addition, we notice that the orbital moments decrease with an increase in the value of U . This is not surprising as the higher values of U localize the Fe 3d electrons more, thereby reducing their hybridization and thus, their orbital moments. Kim et al. [17] had performed X-ray magnetic circular dichroism (XMCD) experiments at the Fe L 23 -edges of GFO and found a net orbital magnetic moment of 0.017µ B at a temperature of 190 K. Assuming that the individual spin and orbital magnetic moments behave as their total magnetization curves, we estimate the net orbital magnetic moment to be 0.034 µ B at 4 K. They estimate their Fe occupancy at the Ga1, Ga2, Fe1 and Fe2 sites to be 0, 0.35, 0.825 and 0.825 respectively. Assuming that the Ga1 and Ga2 moments are parallel to the Fe1 and Fe2 sites respectively, we obtain net orbital moments of 0.007 µ B and 0.006 µ B

for the LDA+U and GGA+U respectively, which are much smaller than experiment. This level of discrepancy is known to be present in LDA/GGA methods and there is no general solution to improve the theoretical orbital moments [53].

B. Origin of Ferrimagnetism in GFO

It is important to point out that the net magnetization obtained theoretically in GFO is zero and hence, we always obtain a perfect antiferromagnetic system. To explain the origin of ferrimagnetism, we performed ab initio calculations to simulate cationic site disorder effects in the smallest possible unit cell. There are three ways that a material with an antiferromagnetic ordering may be ferrimagnetic:

1. The individual magnetic moments are not the same on the two antiferromagnetic sites, thereby giving a net non-zero magnetic moment. Such behavior is seen in magnetite, Fe 3 O 4 , where there are two Fe states, Fe +2 and Fe +3 , which have different moments, and thus give rise to a ferrimagnetic system [54]. In the case of GFO, the Fe ionic state on the cationic sites is Fe +3 .

2. The individual magnetic moments are canted towards a particular direction and give rise to a non-zero magnetic moment along that direction. GFO was thought to be ferrimagnetic due to this reason [12,13], but the experimental observation by Frankel et al. [14] showed that the easy axis of magnetization is along the c-axis.

3. Site disorder in an antiferromagnetic system, where a magnetic atom is replaced with a non-magnetic atom, would leave an uncompensated magnetic moment and the system appears to be ferrimagnetic. Experimental observations on GFO indicate that the Ga and Fe atoms tend to swap places with a preference for the Ga2 site [15]. This might produce an uncompensated moments and hence a ferrimagnetic ground state.

To explain the origin of ferrimagnetism in GFO, we investigated the third possibility concerning site disorder. This is not an easy task because each experimental paper on GFO describes different Fe occupancies at the cationic sites, as can be noticed in the few examples shown in Table II. Moreover, given a set of occupancies, it is not computationally practical to consider larger supercells or even every possible atomic arrangement in the unit cell. However, all experiments show that there is a higher preference of Fe occupying the Ga2 site over the Ga1 site. This could be due to the fact that the Ga1 site is tetrahedral, which would require further investigation to confirm. In addition, the Fe occupancies at the Fe1 and Fe2 sites seem to be identical in all cases. For these reasons, we consider exactly one such case, where the Fe occupancies are as those provided by Hatnean et al.

[33], with Fe@Ga1 = 0 and are possible to simulate in the same cell as earlier. There are obviously many more cells that provide the same occupancies, but those were not considered. Before we present our results, it might be important to first understand the magnetic coupling between the different cationic sites. Given that the Fe1 and Fe2 sites are antiferromagnetically coupled, we have two possible cases for the coupling between the Fe sites and the Ga sites:

1. Case 1: Fe1 is antiferromagnetically coupled to Ga1 and Fe2 is antiferromagnetically coupled to Ga2.

Case 2:

Fe1 is antiferromagnetically coupled to Ga2 and Fe2 is antiferromagnetically coupled to Ga1. A calculation using the occupancies provided in Table II and the magnetic moments for both cases provides us with a net magnetization. Comparing this with the magnetic moments obtained in the corresponding experiments [15,17,33] indicates that the second case is more likely. Using the occupancies as reported by Hatnean et al. and the magnetic moments of 4.5 µ B for Fe, we get the magnetic moments at the Fe1, Fe2 and Ga2 (since the occupation at Ga1 is zero) as -3.375, +3.375 and +2.250 µ B , respectively, as compared to the corresponding experimental values of -4.0, +3.5 and +2.7 µ B . It should be noted that this calculation helps understand the origin of ferrimagnetism in GFO intuitively and does not consider hybridization effects and exchange mechanisms beyond those mentioned, and hence may not be accurate. The disordered unit cell that was used in our ab initio calculation is given in Fig. 3. since they are from different experiments. We have a disagreement of -0.665 µ B and 0.457 µ B for the magnetic moments on the Fe1 and Ga2 sites. However, the magnetic moment on the Fe2 site is in good agreement with experiment as well as the net magnetization, which differs from experiment by -0.619 µ B per unit cell. Moreover, the simple calculation scheme described earlier is a very good indicator of the site magnetization as compared to the ab initio results. These results are similar to those obtained by Roy et al. and thus, we infer that the origin of ferrimagnetism in GFO is due to cationic site disorders.

When we compare the ground state energies between the ideal structure in Fig. 1 and the disordered structure in Fig. 3, we obtain a difference of 231 meV (∼ 2700 K), which is lower than the sum of the predicted values of 25 meV and 400 meV for the Fe2-Ga2 and Fe1-Ga2 disorders respectively as reported by Roy et al. [19]. To study the changes in the electronic structure we plot the total DOS for the ideal and disordered structures in Fig. 4. The figure shows that the two systems are similar in their electronic structure and energy band gaps. Additionally, our calculations indicate that the orbital moments for Fe remain parallel to the corresponding spin moments with a magnitude of 0.017 µ B . Thus, the net orbital magnetization in the system is 0.034 µ B with a direction parallel to the next spin magnetization. Since we considered Fe at the octahedral Fe1, Fe2 and Ga2 sites only, the similarity of DOS and magnitudes of the spin and orbital magnetic moments indicate that the octahedral sites all behave similarly. It might be of interest to study the electronic structure of Fe at the tetrahedral Ga1 site.

Another point that should be made is that the ions were not relaxed and the net forces on the individual atoms increased by a factor of 3 as compared to those in the ideal structure.

Performing ionic relaxations might reduce the ground state energy of the disordered system to the point where the thermal energy (∼ 1000-1500 K) available during synthesis of the experimental samples might be sufficient to cause cationic disorders. These cationic site disorders are statistical and might require larger supercells and more swappings, leading to impractical computational requirements. These are the reasons why ionic relaxations or further investigation of different site disorders were not carried out.

C. Crystal-Field Analysis -Theory, Implementation and Results with the former being lower in energy due to a lower electron repulsion compared to the latter. For a tetrahedral environment, the e g orbitals tend to be lower in energy than the t 2g orbitals due to their relative proximity to the anion p-orbitals. This effect which is purely due to the crystal geometry and can be schematically visualized as shown in Fig. 5.

Crystal

Two important details should be pointed out here, firstly, that the e g and t 2g orbitals are described using real spherical harmonics, also known as cubic harmonics, and secondly, that then the d-orbitals obtained would be a superposition of the maximally split sub-orbitals and can not be clearly defined as e g and t 2g . To obtain the effect of the crystal-field on the e g and t 2g DOS we need either to rotate the real spherical harmonics from the global coordinate system to the local one or make a passive rotation of the coordinate system to align the anions appropriately. We thus have a rotation matrix for the real spherical harmonics analogous to the rotation matrix for the Cartesian coordinate system. Using ab initio techniques, there are at least two ways this can be achieved.

In the first approach, we rotated the crystal structure completely for a cation under study, wherein the wave functions obtained would describe the crystal-field splitting on the given cation appropriately. Such an approach can be used effectively for certain cases only, for example, if there is one cation in a molecular system. In general, however, a crystal does not consist of only one cation surrounded by anions in a given crystal geometry, and thus, the rotations that need to be performed for different cations might be many and different.

In the second approach, we rotated the projected spherical harmonics of the cation under study. This can be done during the same self-consistent cycle and is computationally faster since only the local orbitals are rotated. We have implemented this feature VASP, the details of which are presented in Appendix B). We present below our results for the octahedral cationic Fe1 site in GFO. The crystal-field analysis was performed on all the Fe sites in GFO. Fig. 6 shows the d-projected density of states of an Fe1 atom in GFO, which are similar to all the other Fe sites. The SOC was included in the calculations, which were done using the GGA functional with the Hubbard U = 4 (top) and 8 eV (bottom). The e g and t 2g orbitals are shown in black and red, respectively. We observe that for both values of U , the occupied 3d states are split clearly with the e g orbitals lower in energy than the t 2g . The splitting between the unoccupied states seem to depend upon the value of U , and thus the hybridization,

with the e g orbitals slightly lower in energy than the t 2g for U = 8 eV, and the t 2g orbitals clearly lower in energy for U = 4 eV. Since the Fe1 site (and Fe2 and Ga2 sites as well)

is an octahedral site, the results obtained for the occupied states contradict the splitting expected for octahedral sites. Using the crystal-field splitting, it is possible to probe the filling of d-orbitals and understand why Fe in GFO possesses a non-zero orbital magnetic moment and why the spin magnetic moment of Fe is less than 5 µ B . To do so, we plot the ( , m, s)-projected DOS of the 3d electrons of the same atom in Fig. 7. Also shown are the effects of the SOC on the DOS. This is an Fe1 site which is why the majority spins are of spin-down nature.

Moreover, there are some spin-up states, mainly of t 2g character, which belong to the neighboring oxygen atoms. These states oppose the moments on the Fe sites and create a small non-zero orbital moment and reduce the spin magnetic moment from its maximum value of 5 µ B . This would indicate that the bonding between the Fe and O ions in GFO has a partial covalent character. As a result of this covalent character, certain spin-polarized states should be present on the oxygen ions. Our calculations indeed agree with this prediction, indicating small spin magnetic moments between 0.010 µ B and 0.087 µ B for U = 4 eV, and between 0.008 µ B and 0.065 µ B for U = 8 eV on the oxygen atoms. These values are more affected by the level of hybridization controlled by U than the SOC. We thus conclude that increasing the value of U drives the Fe system to a +3 state and that the bonding between the Fe and O atoms is partially ionic. Furthermore, our results show that this bonding is not influenced by the SOC, with the exception of creating a non-zero orbital moment. The conclusions drawn here support the calculations done by Ibrahim and Alouani [18], where they obtained Bader charges of +1.63 to +1.73 on the Fe atoms depending on U , instead of +2 as expected for an ionic crystal, and concluded a partial ionic bonding.

To understand the mechanism causing a tetrahedral-type splitting in an octahedral environment, a point-charge model was developed, in which the effects of hybridization can be fully neglected. Such a model would help us understand the e g -t 2g splitting obtained for the two different values of U . Qualitatively, the electrons in the model can only interact through the Coulomb interaction with the neighboring negatively charged oxygen ions. With this requirement in mind, we assumed ionic charges of +3 on the central atom and -2 on the neighboring atoms forming the octahedron. For maintaining consistency between the ab initio calculation and the point-charge model, the octahedron was oriented in the same way as in the rotated frame of reference, i.e. with the octahedral arms aligned maximally along the axes. We then rewrite the Coulomb potential felt at the central atom in terms of the spherical harmonics as [START_REF] Jackson | Classical Electrodynamics[END_REF],

V (r) = 6 i=1 q i |r -R i | = 4π 6 i=1 q i ∞ =0 m=- (-1) m 1 2 + 1 Y -m (r) Y m ( Ri ) r < r +1 > , (1)
where q i = Ce 2 is an effective charge of electron-electron interaction, which describe the interaction of the ligand effective charge of Ce with an electron of charge e in the 3d. Here we set the effective charge on the neighboring ligand atoms to two, which corresponds to that of O 2-in GFO. The Y m is the spherical harmonics, and R i is the distance vector connecting the central atom to the ligand i. Note that the local reference frame is centered on the iron atom and its axes are along the ligands if the octahedron is perfect. Since the octahedron is deformed, the local frame is optimized to reduce the angles between R i and the corresponding axis (see appendix A for details). The spherical angles θ and ϕ of the vector R i are then obtained in this optimized reference frame to compute the cubic harmonics Y m ( Ri ). Using the above expression, we can calculate the matrix elements M m,m of the ligand's potential between the central iron 3d orbitals as,

M m,m = φ =2 Y m =2 |V |φ =2 Y m =2 = 4π 6 i=1 q i ,m (-1) m+m 2 + 1 R i Y m ( Ri ) Y -m 2 (r) Y -m (r) r < r +1 > (φ 2 (r)) 2 Y m 2 (r) dr, (2)
where φ 2 (r) is the radial function corresponding to the Fe 3d orbitals. The above expression contains a product of three spherical harmonics and their integration, which can be simplified using the Gaunt coefficients, which in terms of the Wigner 3-j symbols are,

M m,m =5 6 i=1 q i ,m (-1) m+m 4π 2 + 1   2 2 0 0 0     2 2 m -m -m   Y m ( Ri ) ×   1 R +1 i R i 0 r +2 (φ 2 (r)) 2 dr + R i ∞ R i (φ 2 (r)) 2 r -1 dr   , (3)
where

  j 1 j 2 j 3 m 1 m 2 m 3 
 are the Wigner 3-j symbols linked to the Clebsch-Gordon coefficients as,  

j 1 j 2 j 3 m 1 m 2 m 3   = (-1) j 1 -j 2 -m 3 √ 2j 3 + 1 j 3 , -m 3 |j 1 , m 1 ; j 2 , m 2 . ( 4
)
The Wigner 3-j symbols are used for their symmetry properties, which allow us to restrict the summation over , m in Eq. ( 3) to = 0, 2, 4 and m = m -m.

Eq. ( 3) can be computed directly using the atomic 3d wave function φ 2 of Fe as obtained by an all electron atomic program, however in our case since the crystal field is an empirical method used here only to understand our ab initio calculations, we rewrote the matrix elements in terms of the hydrogen-like radial function for the 3d orbitals, corresponding to the 3d orbitals of Fe. Instead of using the atomic number Z = 26 we used an effective Z = 6.25 as suggested by the so called Slater's rules [START_REF] Tudela | [END_REF]. According to Slater, electrons within the same group of 3d electrons shield 0.35 charge, whereas electrons with lower orbital groups shield 1 charge. For iron the total shielding experienced by a 3d electron is 19.75 charges, which led to Z = 6.25. The electrons in the higher 4s orbital do not contribute to the shielding of the nucleus of iron. The use of the hydrogenoid wave function led to the following matrix elements:

M m,m = Z 216a 0 6 i=1 q i =0,2,4 m (-1) m+m Y m ( Ri ) 4π 2 + 1   2 2 0 0 0     2 2 m -m -m   , × 1 S +1 i γ(l + 7, S i ) + S i Γ(6 -l, S i ) (5)
where S i = 2Z R i 3a 0 , and where a 0 = 4πε 0 me 2 is the Bohr radius. Here γ(l, x) and Γ(l, x) are, respectively, the lower and upper incomplete gamma functions. They can be calculated recursively as γ(l, x) = -x -1 e -x + (l -1)γ(l -1, x) and Γ(l, x) = x -1 e -x + (l -1)Γ(l -1, x)

starting from γ(1, x) = 1 -e -x and Γ(1, x) = e -x ,
For the 3d orbitals, we have five m values and thus, M is a 5 × 5 matrix. Note that the spin degrees of freedom are not included. This is justified by the fact that Fe in GFO has essentially a +3 ionic charge and a large band gap, where the states of one spin interact mainly among themselves than with the other spin states. However, if one wishes to include the spin quantum number and the SOC, it can be done using a 10 × 10 matrix, where the off-diagonal 5 × 5 matrices would contain the couplings (exchange, SOC) between the two spin states.

Before diagonalizing the matrix M , we made a unitary transformation from the spherical harmonics to cubic harmonics, so that the ensuing eigenvectors are directly expressed as a linear combination of e g and t 2g states. We set the order of the cubic harmonics basis as t 2g (xy, yz, and zx) and then e g (3z 2 -r 2 , x 2 -y 2 ). The new matrix to diagonalize U † M U is directly in the cubic harmonics basis set, where U is the unitary transformation between the cubic and spherical harmonics.

The results were calculated for the same Fe1 atom as described in Fig. 6, and are plotted in Fig. 8 together with the weight of each eigenvector. The figure shows that there is very low mixing of the e g and t 2g orbitals. We obtain a similar behavior as obtained for the unoccupied spin-up states using the ab initio splittings, with the t 2g orbitals lower in energy than the e g orbitals, but not for the occupied spin-down states. Since the pointcharge model is essentially based on unscreened Coulomb interaction with point-like ionic neighbors without hybridization, the crystal-field splitting of each t 2g or e g orbital arises only from the distorted octahedral environment. There are at least two inferences we can draw from the results obtained in this section,

• The conventional meaning of e g and t 2g orbitals is valid because the crystal-field splitting is in agreement with what is expected (Fig. 5). However, even though it reproduces the splitting of the spin-up DOSs, it does not agree with e g and t 2g splitting of the occupied spin-down of the Fe1 atom (see Fig. 6).

• Although the octahedral environment is so distorted, it resulted only on slight splitting of the two e g and three t 2g as seen in Fig. 8. However the crystal-field splitting between the e g and t 2g of 1.6 eV is in agreement with the ab initio results for the empty states of 1 eV (for U = 4 eV) but less with the DOSs for U= 8 eV.

While the first point seems right, it is possible to distinguish clearly the e g and t 2g orbitals of the unoccupied states of the ab initio calculation of 1 eV (for U = 4 eV) as well as the point-charge model 1.6 eV, due to the low mixing of the orbitals. However for the occupied spin-down states, the point charge model failed completely, as it can not explain how the e g states become lower than the t 2g . One can not just assume that the spin-up and spin-down DOSs are rigidly split by 8 eV (for U= 4 eV) and 11 eV (for U= 8 eV) by the exchange interaction. The hybridization between the oxygen p states and the iron e g and t 2g is key for understanding the physics of GFO beyond the point-change model.

As indicated above, the second point is less intuitive. Even though the Fe central atom is displaced away from the center causing a high amount of distortion of the octahedral environment, as was first shown in an earlier work [18], it resulted only in small splitting and a small mixing of the e g and t 2g orbitals.

The point-charge model is not an exact model and can only provide a qualitative understanding of the crystal-field splitting for systems which are mainly ionic and the hybridization between the central atom orbitals and neighboring ligands is low. As the amount of hybridization increases (lowering of U ), the model is less appropriate to study the orbital splitting behavior, since the neighboring ions can no longer be considered as fully ionic or point-charge like. However, the point-charge model remains useful and is often used for estimating the splitting of the d or f orbitals.

IV. MAGNETOELECTRIC PROPERTIES

GFO has been known as a ferrimagnetic and ferroelectric material, i.e. it contains a spontaneous magnetization and a spontaneous electric polarization. Additionally, these two ferroic orders have been known to be coupled to one another [8, 11-13, 15, 58]. However, a theoretical understanding of the mechanism driving this coupling is missing and until recently, even the electric polarization of the GFO system was not computed using ab initio techniques [31,32]. With an aim of understanding the coupling between the magnetic and electric ferroic orders, we performed ab initio calculations for the electric polarization and its dependence on the magnetization direction as implemented in VASP [46][47][48][49][START_REF] Resta | Berry Phase in Electronic Wavefunctions[END_REF][51].

A. Temperature Dependence of Electric Polarization

To simulate a temperature dependence, we assume a linear interpolation of the atomic positions and lattice vectors between those at 4 K and 230 K. The initial (at 4 K) and final IV. It can be observed from this data that there is very little change in the atomic positions and the lattice vectors, for example, change in volume is less than 0.2 %, and thus, we can work in the linear regime between these temperatures.

To simplify work in the linear regime, we defined a parameter λ between 0 and 1 defining the temperature and positions at 4 K and 230 K respectively. For example, the temperature is then defined as,

T (λ) = 4 + (230 -4)λ (6)
λ was varied in steps of 0.1 and we calculated the electric polarization for the atomic positions at these points. Since the initial and final positions are obtained from experiment, the ionic relaxation was not carried out. Additionally, the electric polarization is multi-valued differing by integral values of the polarization quanta [49,59]. To resolve this, it is necessary to calculate the polarization on a path connecting the non-polar structure (centrosymmetric) and the polar structure (non-centrosymmetric). As a result, the electric polarization values were already manipulated based on the method developed by Stoeffler [32], who had performed the calculation for the GFO structure at 4K along a polarization lattice branch connecting the corresponding non-polar and polar structures.

The absolute values of the electric polarization for the ideal system as a function of temperature are shown in Fig. 9. The polarization vector is aligned along the negative y-axis, with a magnitude of about 23.5 µC/cm 2 , close to the value of 25 µC/cm 2 as reported by Stoeffler [32]. It is seen that as the temperature increases, the magnitude of the polarization decreases. This is expected since increase in temperature leads to an increase of the inter-atomic distances, which causes hybridization to decrease and electrons to be more localized near the parent atom and the bonding becomes less ionic. Since the electric polarization is a measure of how far apart the charges are, a localization of electrons near the ions implies a lowering of the magnitude of the electric polarization.

This can be understood with the help of Fig. 10, which shows the electronic and ionic components of the electric polarization in GFO. The polarization is plotted along the posi- tive y-axis, which is why the net polarization is negative. The ionic part of the polarization increases in magnitude with temperature. Since, the ionic part of the polarization is the sum of the dipole moments per unit cell, an increase of temperature would drive the atoms apart, thereby increasing the dipole moments and thus, the electric polarization.

The electronic part of the polarization also increases with temperature, however in the opposite direction, due to the sign of the charge. If we consider that the Wannier centers of the electrons are located at the positions of the ions, then increasing the temperature would drive the ions apart and the Wannier centers, hence increasing the polarization.

However, if this were the only factor, we would expect no change in the net polarization since the ionic positions and the Wannier centers would be driven apart by exactly the same amount. The other factor is the electron hybridization, which decreases with increase in the inter-atomic distances. As a result of this, the electronic density becomes more localized and the bonding becomes more ionic, which increases the electric polarization.

Since the electronic polarization and the ionic polarization are in opposition, the net polar- The electric polarization obtained in the previous subsection are in good agreement with what was found by earlier theoretical studies [32]. Recent experimental studies [60] were performed on Ga 2-x Fe x O 3 with x = 1.1, which indicate polarization values around 33 µC/cm 2 , not very far from ours. These confirmations support our calculations and allow us to study magnetoelectric effects in GFO. Magnetoelectric effects form a class of phenomena that arise due to the coupling between the magnetic and electric ferroic orders, which are the ferrimagnetic and ferroelectric orders in GFO. In our work, we probe the effect of rotation of the magnetization axis on the electric polarization.

There are two types of magnetoelectric effects -direct and indirect. The direct magnetoelectric effect arises from the SOC, that couples the spin and the lattice. When an external magnetic field is applied, the electrons move to a different ground state that causes the electronic polarization to change. This effect does not require the ionic positions to change and is a direct consequence of the external magnetic field, hence the term direct. On the other hand, the indirect magnetoelectric effect is a consequence of the external magnetic field moving the ions and altering the volume of the cell, thus changing the electric polarization.

This effect is also called the magnetostrictive magnetoelectric effect. Based on experiments conducted in the 1960's, it was hypothesized that the magnetoelectric effects observed in GFO are due to the indirect mechanism [12]. However, in the 1990's, Popov et al. indicated that the direct mechanism is responsible for the observed magnetoelectric effects in GFO [61]. In spite of the unresolved problem for over half a century, there has been no systematic theoretical study trying to probe this phenomenon to the best of our knowledge.

We begin by rotating the magnetization axis in the y-z plane, starting from the +z-axis, the original configuration, to the -z-axis through the +y-axis. Rotation of the magnetization axis is akin to rotation of the moments such that they remain parallel (or anti-parallel)

to the magnetization axis. Since we can not include external magnetic fields in our calculation, we fix the moments along a rotated magnetization axis to simulate a saturation external magnetic field in the y-z plane that will rotate the moments in the same way. By maintaining an antiferromagnetic system, we assume that the exchange energy, i.e. the energy to overcome the antiferromagnetic coupling between the Fe ions, is always higher than the energy of the corresponding external magnetic field. To check whether the direct mechanism is responsible for the observed magnetoelectric effects, the calculations were performed keeping the ions fixed. As a result, any change in the electric polarization is due to the change in the electronic polarization only. While a clear pattern exists in both curves, the change is less than 0.4 nC/cm 2 , five orders of magnitude lower than the spontaneous electric polarization determined in the previous subsection at 23.5 µC/cm 2 . Moreover, it is very hard to measure these small values experimentally. Thus, quantitatively speaking, the direct mechanism is far from sufficient to explain the magnetoelectric effects in GFO. The indirect mechanism on the other hand might be the main cause for the observed magnetoelectric effects. This is supported by the temperature dependence calculation results, where the volume change of about 0.2 % caused the electric polarization to change by about 0.6 µC/cm 2 or 2.5 %.

To properly determine the main mechanism for the magnetoelectric effects, it is hence important to perform similar calculations while allowing the unit cell volume and ionic positions to change, which might be very expensive since GFO is a relatively large system.

These calculations can be done by performing only a volume relaxation to minimize the computational effort, and have not yet been carried out to the best of our knowledge.

V. CALCULATION OF XAS AND XMCD

X-ray absorption spectroscopy (XAS) and x-ray magnetic circular dichroism (XMCD) are excellent methods for probing the electronic and magnetic structure of materials. These methods have been under intense study, both theoretically and experimentally [START_REF] Ebert | Electronic Structure and Physical Properties of Solids[END_REF][63][64][65][66][67][68][69][START_REF] Van Der Laan | Magnetism and Synchrotron Radiation: Towards the Fourth Generation Light Sources[END_REF][71]. With XAS and XMCD, it is possible to obtain element-specific structure and magnetization. The chemical selectivity of each core orbital of any atomic species in a material is what makes XAS and XMCD more capable of characterizing magnetic systems than traditional magnetic techniques. This property is used to study magnetism in 3D, 2D and 1D systems. XMCD originates from the coupling between the photon helicity (± ) and the atomic magnetic moments, thereby creating a difference between the absorption cross sections measured with respect to the magnetization axis for left and right circularly polarized light. Under the PAW scheme, the x-ray absorption σ at frequency ω for a given polarization µ is given by

[72] σ µ (ω) = 4πα 3 m 2 ω M,n,k,s p, ,m,m C J,M ,m ;1/2,s C ,m 1,µ; ,m C ,0 1,0; ,0 , 0|∇ 0 |p, , 0 P n,k,s p, ,m 2 δ( ω -nks + JM ), (7)
where, p is the projector index, , m are the angular momentum and magnetic quantum numbers for the valence states, , m are those for the core states, J, M is the core spinorbit-split quantum numbers and C J 1 ,M 1 J 2 ,M 2 ;J 3 ,M 3 = J 1 , M 1 |J 2 , M 2 ; J 3 , M 3 are the Clebsch-Gordon coefficients. P n,k,s p, ,m are the projectors in the PAW method represented in the complex spherical harmonics basis (see Ref. [72] for details of our implementation).

The calculated L 2,3 XAS edges and XMCD for GFO with disorder and for U = 8 EV are provided in Fig. 12. The theoretical and experimental spectra for the individual circularly left and right polarized light were not discernible, indicating that the system is essentially antiferromagnetic. Indeed we observe that the two curves are very close to each other.

The computed XAS was then averaged and then normalized to match the experimental L 3 highest peak, and the L 2 -peak was shifted by the 13 eV spin-orbit splitting of the 2p 1/2 and 2p3/2 core states in agreement with the position of the experimental counterpart. Fig. 12 shows that the experimental peaks exhibit multiplet splitting. This behavior occurs generally in 3d and 4f oxides and compounds because the 3d or 4f electrons are strongly correlated and localized. DFT is a single-particle picture and fails to correctly account for the splitting of localized orbitals, but overall the calculated DFT spectra are in good agreement with experiment. The multiplet splitting can be described using atomic multiplet models [73].

The XMCD signal was also computed after normalizing the XAS and shows a dichroism indicating a net magnetic moment in the system. The sum rules [67,68] were used to compute the average orbital and spin moments per Fe atom in the system and are shown the XMCD sum rules compared to the experimental results [17] in Table V. The number of holes was set to 5, since Fe in GFO has an official ionic charge of +3, which amounts to a half-filled 3d shell. The computed spin moment agrees well with experiment, while the orbital moment is underestimated by a factor of 2. Note however that the sum rules should in principle provide the same average magnetic moment per iron atom as the direct ab initio calculation. However, the ab initio calculation for U = 4 and 8 eV (see table III) show that the spin magnetic moment per iron atom are respectively 1.06 and 1.13 µ B , and the average orbital moment of 0.006 and 0.0055 µ B . Thus the sum rules produced an error of -38% for the spin moment and +45% increase for the orbital moment. These types of error due to the use of the XMCD sum rules are well documented in the literature [69].

Additional source of errors could also be due to the fact that the PAW 2p core state is frozen and the basis set includes only the augmentation region contribution. If the experimental errors form the sum rules are the same as for the calculated values then the experimental spin moment per atom should be 1.4 µ B in agreement with the experimental averaged moment per iron atom of 1.2 µ B as produced from table III. As for the average experimental orbital moment, taken into account the sum rule errors, it is expected to be 0.011 µ B , which is twice the theoretical value. Note that, as indicated above, DFT underestimates the orbital moment by a factor of two [53]. In conclusion, the overall agreement with experiment is satisfactorily, since without any site disorder GFO will not exhibit any XMCD signal.

VI. DISCUSSION

In this work, we presented the results of our ab initio calculations on GFO. To correctly account for the strongly correlated 3d electrons the LDA+U model was used, where we found that the value of U = 8 eV best reproduced the experimental values of the energy band gap and the site magnetizations. We also found that the inclusion of the SOC did not change the electronic properties drastically and allowed us to obtain a non-zero orbital moment in GFO.

By performing site disorder studies based on cationic occupancies we confirmed earlier hypotheses attributing these disorders to the cause of ferrimagnetism in the system. The difference between the ground state energies of the ideal and disordered systems is below what was previously estimated and much closer to the thermal energy available during synthesis of the system. This indicates that site disorder can not be completely controlled using current synthesis methods. Further studies using excess-Fe might provide clues to better control the electronic and magnetic properties, in spite of the site disorders, to help the development of practical devices.

To understand the magnetic ordering of the Fe sites, we implemented the crystal-field analysis in VASP. For the octahedral Fe sites, we obtained the t 2g states below the e g states as expected from an octahedral splitting. However, this model calculation was found to be in contradiction with our ab initio results as shown in Fig. 6 for the occupied states. This

shows that the ab initio interactions are more complex than these given by a point charge model. As stated in crystal-field section, the hybridization between the oxygen p states and the iron e g and t 2g states is key for understanding the physics of GFO beyond the pointchange model. Indeed when the extended oxygen p state inter the Fe augmentation region it is decomposed in the local spherical harmonics and looks like a d state. This makes the interpretation of the Fe d states very difficult as it is partly composed of oxygen states. We also showed that the chemical bonding between Fe and oxygen is partly covalent, thereby explaining the non-zero orbital moment and supporting earlier Bader charge calculations.

This doesn't say that iron is not in Fe 3+ , it only claims that the oxygen O 2-Wannier function is very extended, i.e., the two electrons gained by the oxygen are more spread contrary to the assumed formal charge picture like, for example, in NaCl.

In addition, we presented our results for the linear temperature dependence of the electric polarization. The polarization has a monotonic decay which could be explained with the simple model of reduction of hybridization due to increase of inter-atomic distances. This was interpreted by separating the electronic and ionic parts of the electric polarization.

We also presented the results of our direct magnetoelectric studies, where we rotated the magnetization axis to simulate a rotation of the moments due to an external field while keeping the ions and cell volume fixed. The maximum change in the electric polarization was under 0.4 nC/cm 2 , and indicates that the direct mechanism fails to explain the observed magnetoelectric effects in GFO. To truly find the cause of these effects, one must perform similar studies of the electric polarization by allowing the cell volume to change due to the change of the magnetization axis, i.e. a magnetostriction induced polarization.

Finally, we presented the results for the calculation of the XAS and XMCD for the Fe L 2,3 -edges in GFO. The theoretical and experimental spectra match, indicating that most of the features are captured by the calculations. The lack of multiplets in the theoretical curves indicate the shortcoming of the single-particle picture, and further work needs to be done to capture them using ab initio techniques. The sum rules for the spin moments are in good agreement, and similar studies can be carried out for GFO with different site disorders for comparing theoretical calculations with experiment.

VII. CONCLUSIONS AND PERSPECTIVES

In summary, the results for the ab initio calculations on GFO performed under the LDA+U and GGA+U approximations with SOC were presented. The value of U =8 eV best reproduced the experimentally observed values of the energy band gap and the site magnetizations. U was shown to push the 3d electrons away from the Fermi level and thus played an important role in the electron hybridization. The inclusion of SOC did not affect the electronic structure of the system. Cationic site disorder studies were then performed using experimental cationic occupancies, which reproduced the experimentally observed ferrimagnetism. The energy difference was shown to be closer to the thermal energy available during synthesis, thereby indicating difficulty of site disorder control.

We have implemented the crystal-field analysis in VASP and used it to show that the tetrahedral-like splitting of the Fe occupied e g and t 2g DOS obtained in the ab initio calculation can not be explained by a point-charge interaction. Indeed, the interaction between the oxygen ligands and the Fe d states is very complex and it beyond the simple point-charge model, as the hybridization plays a key role in the e g and t 2g splitting. This is in agreement with the partly covalent bonding between Fe and O, and explain the origin of the non-zero orbital moment.

We then presented the linear temperature dependence of the electric polarization, which has a monotonic decay explained as separation of the ions with temperature. Then, we rotated the magnetization axis of the system to simulate a rotation of the magnetic moments under an external field. By fixing the ions and cell volume, the direct magnetoelectric effect was computed. The values obtained were very small indicating that the indirect mechanism (magnetostrictive induced polarization) might be the main cause for the observed magnetoelectric effects in GFO.

Finally, we calculated the XAS and XMCD for the Fe L 2,3 -edges in GFO with site disorders and compared them with experiment. We found good agreement with experiment for the spectra as well as the sum rules.

As perspectives for the system, it would be fruitful to study excess-Fe along with cationic site disorders to better control the electronic and magnetic properties of GFO. The crystalfield analysis can be redone on the Ga sites, both octahedral and tetrahedral, to understand the magnetic ordering of disordered systems. Magnetoelectric effects of the indirect nature can be studied to help understand the coupling between the ferroic orders. The methods developed here can also be applied to other multiferroic systems, such as Cr 2 O 3 or BiFeO 3 .

VIII. ACKNOWLEDGMENTS

This work was performed using HPC resources from the Strasbourg Mesocenter and from the GENCI-CINES Grant gem1100. We thank Nathalie Viart and Bohdan Kundys for a useful discussion on magneto-electric coupling from an experimental perspective.

We define the corresponding 5 × 5 rotation matrix  for the real spherical harmonics with 2. m = -1: ).

a 21 =
We now derive the form of A for rotations about each of the Cartesian axes. Since any rotation can always be represented in terms of rotations about the Cartesian axes, R z , R y and R x , the A matrices can be used for all cases. As is shown, the double solution obtained for a i3 are always identical and thus, all rotations of the spherical harmonics can also be represented in terms of A z , A y and A x .

We begin by considering a rotation about the z-axis. In the real spherical harmonics basis, Rz (θ) → A z (θ). The R matrix in Cartesian basis is,

R z (θ) =      cos θ -sin θ 0 sin θ cos θ 0 0 0 1      . ( B3
)
It can be easily verified that the double solutions for a i3 are in fact identical for this case (and for the following two cases as well). Thus, the A z (θ) matrix is, A Z-Y -Z (α, β, γ) = A z (γ)A y (β)A z (α). The order is important since matrix multiplication is not commutative. In our implementation, we use the A Z-Y -Z (α, β, γ) matrix directly after computing the Euler angles α, β and γ of each octahedron.

A z (θ) =          

FIG. 1 .

 1 FIG. 1. (Color online) GaFeO 3 unit cell with Fe, Ga and O atoms in yellow, red and blue respectively.

FIG. 2 .

 2 FIG. 2. (Color online) GGA+U and LDA+U calculated total DOS for GFO for U = 4 and 8 eV.

Property

FIG. 3 .

 3 FIG. 3. (Color online) Unit cell simulating a disordered configuration to understand the origin of ferrimagnetism in GFO.

  [19], who predicted that the individual Fe1 and Fe2 moments do not change by a large amount with respect to the perfect structure. However, they observed that the Fe ion at the Ga2 site has a relatively lower moment of 4.11 µ B , in contrast to our observed moment of 4.486 µ B . The corresponding individual moments of the Fe ion at the Ga2 sites based on experimental average values of 2.7 µ B are 5.4 µ B , which does not seem likely considering that the maximum magnetic moment Fe can possess is 5 µ B . The same is true for the value of -4.0 µ B reported for the Fe1 site, which corresponds to individual moments of -5.33 µ B based on the iron occupancy of 0.75. The net magnetic moment obtained is about 10 µ B per unit cell,

FIG. 4 .

 4 FIG. 4. (Color online) Comparison of the total DOS for the ideal and disordered structures of GFO for U= 8 eV. Similar results are obtained for U=4 eV.

  -field theory helps describe the splitting of the electron orbitals of an atom, usually a d or f cation, in the presence of a Coulomb potential generated by neighboring atoms, usually anions. This model has been very successful in analyzing d and f splitting for different magnetic orders and when combined with molecular orbitals, has successfully explained spin crossover phenomena, where the energy gap of the splitting can be of the order of the pairing energy between the electrons[55]. In GFO we have two types of crystalfield environments in the system -octahedral and tetrahedral. An Fe ion at the center of these environments would be influenced by six and four O anions, respectively, and the splitting of the degenerate 3d orbitals of Fe would reflect this interaction. In the case of an octahedral environment, the electrons of the anions are closer to the the d z 2 and d x 2 -y 2 orbitals of Fe and away from the d xy , d yz and d xz orbitals. The strong electron-electron repulsion would lead to a splitting of the five d-orbitals into two sub-shells, called t 2g , consisting of the d xy , d yz and d xz orbitals, and e g consisting of the d z 2 and d x 2 -y 2 orbitals,

FIG. 5 .

 5 FIG. 5. (Color online) Crystal-field splitting of 3d orbitals in octahedral and tetrahedral fields

1 FIG. 6 .

 16 FIG.6. (Color online) Crystal-field effect on the Fe1 e g (black) and t 2g (red) DOS using GGA+U method for U= 4 eV (upper panel) and 8 eV (lower panel).

FIG. 7 .

 7 FIG. 7. (Color online) Effect of SOC and U on the the spin polarized (positive values spin up and negative values spin down) Fe1 e g (black) and t 2g (red) DOS

FIG. 8 .

 8 FIG. 8. (Color online) Point-charge model crystal-field splitting of the 3d states of iron in the FeO 6 octahedron environment in GFO into e g and t 2g states. The different colors represent the weight in percent for each of the basis functions of the e g and t 2g symmetries. The total weight for each state adds to 100%.

(

  at 230 K) positions were measured with neutron diffraction patterns by Arima et al. and are given in Table

K with a = 8 .

 8 71932 Å, b = 9.36838 Å and c = 5.06723 Å at 4 K and a = 8.72569 Å, b = 9.37209Å and c = 5.07082 Å at 230 K.[15]

FIG. 9 .FIG. 10 .

 910 FIG. 9. (Color online) Calculated electric polarization in GFO as a function of temperature. The dashed line is a linear fit.

6 ∆P (nC/cm 2 )

 62 FIG. 11. (Color online) Change in the electric polarizations along the y-and z-axes in GFO as a function of rotation of the magnetization axis

FIG. 12 .

 12 FIG. 12. (Color online) Calculated left and right circularly polarized L 2,3 -edges at the Fe1 site (left) and averaged (right) XAS and XMCD of Fe in disodred GFO, calculated for U = 8 eV, compared to experiment[17].

= 2 ,aa 11 = 21 a 12 = 21 a 15 =

 21121122115 and by manipulating the indexes a ij , we can write Y m (r ) = m a mm Y m (r). If we write out the full form of the real spherical harmonics with r i = j R ij r j and use |r | = |r|, are required and not the final positions themselves. However, the final positions help us understand the amount of distortion in the geometry, and are important while developing the point-charge model, where the Coulomb potential due to the anions is expanded as a series of spherical harmonics. The all-electron wave function in the projector augmented wave basis (PAW) is given as, |n, k, s = |n, k, s + p, ,m |p, , m -|p, , m P nks p m , (A3) where |n, k, s is the pseudo-wave function for band n, wave vector k and spin s. |p, , m and |p, , m are the partial and pseudo-partial waves corresponding to projector type p, angular momentum quantum number and magnetic quantum number m. P nk p m are the projector coefficients. In this equation, SOC is not included, which is why the spin quantum number s appears explicitly. The inclusion of SOC is trivial and does not affect the final results in our work. The local site projected DOS is defined within the augmentation region, where the all-electron wave function is given only by the partial waves as, |n, k, s = p, ,m,s P nks p m |p, , m . (A4) When the Cartesian reference frame is rotated, the new real spherical harmonics maintain their quantum number and the total DOS for a given remains unchanged. Hence, we have, |Ψ nk = p,m,s P nks p m |p, , m = p,m,s P nks p m |p, , m , (A5) where the 'barred' objects are the rotated quantities. To find how the new projector coefficients transform, we express the rotated real spherical harmonics in terms of the original basis set. This is trivial since the radial functions do not change, and only the spherical harmonics are rotated. Thus, |p, , m = m a mm |p, , m . (A6) When we substitute this in Eq. (A5), we get, mm |p, , m = p,m ,s m a mm P nks p m |p, , m = p,m ,s P nks p m |p, , m . R 11 R 22 + R 12 R R 12 R 23 + R 13 R 22 a 13 = √ 3R 13 R 23 = -√ 3(R 11 R 21 + R 12 R 22 ) a 14 = R 11 R 23 + R 13 R R 11 R 21 -R 12 R 22 .

  a general rotation about the y-axis. In the real spherical harmonics basis, Ry (θ) → A y (θ). The R matrix in Cartesian basis is, the R matrices, the A matrices can be multiplied in the order of rotations. This is a faster way to generate the net rotation matrix since the individual matrices are simpler. One can of course, first create a net R matrix by multiplying the individual R matrices and then using the A matrix element definitions to create the A matrix, i.e.

  Table III shows our average site magnetization results with ab initio techniques and compares them to the experimental and the simple calculation results.

	m S @Fe1	-4.0	-3.01	-3.09	-3.34	-3.38
	m orb @Fe1		-0.017 -0.017	-0.012 -0.013
	m S @Fe2	3.5	3.08	3.09	3.35	3.38
	m orb @Fe2		0.018 0.017	0.014 0.013
	m S @Ga1	0.0	0.0	0.00	0.00	0.00
	m orb @Ga1	0.0	0.00	0.00	0.00	0.00
	m S @Ga2	2.7	2.05	2.06	2.24	2.26
	m orb @Ga2		0.011 0.011	0.009 0.008

The table shows that the ab initio super cell site disorder calculations agree relatively well with experiment. Note that the experimental values are different from those in table I U = 4 eV U = 8 eV m S /m orb (µ B ) Expt.

[33]

Disorder Average Disorder Average TABLE III. Ab initio spin (m S ) and orbital (m orb ) averaged magnetic moments per site for a disordered GFO (Disorder) compared with experiment and simple calculation based on site occupations (Average) as given in the middle column of table II.

TABLE IV .

 IV Experimental positions in units of the lattice parameters of GFO ions at 4 K and 230

	Site	4 K			230 K
	x	y	z	x	y	z
	Ga1 0.1500 0.0000 0.1781 0.1501 0.0000 0.1761
	Ga2 0.1593 0.3073 0.8106 0.1597 0.3067 0.8091
	Fe1 0.1538 0.5831 0.1886 0.1525 0.5827 0.1893
	Fe2 0.0346 0.7998 0.6795 0.0351 0.7992 0.6787
	O1 0.3228 0.4262 0.9716 0.3223 0.4260 0.9740
	O2 0.4864 0.4311 0.5142 0.4877 0.4313 0.5168
	O3 0.9979 0.2022 0.6541 0.9963 0.2008 0.6521
	O4 0.1593 0.1974 0.1480 0.1593 0.1961 0.1475
	O5 0.1695 0.6717 0.8437 0.1715 0.6714 0.8410
	O6 0.1736 0.9383 0.5166 0.1725 0.9379 0.5153

TABLE V .

 V Calculated iron-site averaged orbital m orb and spin m s moments (in units of µ B ) from

	Magnetic Moment Theory Experiment [17]
	m s	0.777	0.870
	m orb	0.008	0.017

  R 21 R 32 + R 22 R 31 a 22 = R 22 R 33 + R 23 R 32 a 23 = √ 3R 23 R 33 = -√ 3(R 21 R 31 + R 22 R 32 ) a 24 = R 23 R 31 + R 21 R 33 a 25 = R 21 R 31 -R 22 R 32 . R 32 -R 11 R 12 -R 21 R 22 ) R 33 -R 12 R 13 -R 22 R 23 ) R 33 -R 11 R 13 -R 21 R 23 ) a 41 = R 31 R 12 + R 32 R 11 a 42 = R 32 R 13 + R 33 R 12 a 43 = √ 3R 33 R 13 = -√ 3(R 31 R 11 + R 32 R 12 ) a 44 = R 33 R 11 + R 31 R 13 a 45 = R 31 R 11 -R 32 R 12 .5. m = 2:a 51 = R 11 R 12 -R 22 R 22 a 52 = R 12 R 13 -R 23 R 23 R 13 R 11 -R 23 R 21

	a 53 =	√ 2 3	(R 2 13 -R 2 23 ) = -	√ 2 3	(R 2 11 -R 2 21 + R 2 12 -R 2 22 )
	a 54 = a 55 =	1 2	(R 2 11 -R 2 21 -R 2 12 + R 2 22
	3. m = 0:
	a 31 = (2R 31 a 32 = 1 √ 3 1 √ 3 (2R 32 a 33 = 1 2 (2R 2 33 -R 2 13 -R 2 23 ) = -	1 2	(2R 2 32 -R 2 12 -R 2 22 + 2R 2 31 -R 2 11 -R 2 21 )
	a 34 = (2R 31 a 35 = 1 √ 3 1 √ 3 (2R 2 31 -R 2 11 -R 2 21 -2R 2 32 + R 2 12 + R 2 22 ).
	4. m = 1:

Rotation of the real spherical harmonics can be done in two ways -firstly by rotating the spherical harmonics using the Wigner D-matrices and then rewriting them as the real spherical harmonics, and secondly by rotating the real spherical harmonics directly by using the rotation matrix of the coordinate system. The Wigner D-matrices are well described by Bradley and Cracknell and the interested reader is directed to their work [START_REF] Bradley | The Mathematical Theory of Symmetry in Solids[END_REF]. In our implementation, the second method is used.

Our implementation involves rotation of the d-orbitals following a rotation R(α, β, γ) of the Cartesian coordinate system. Thus, we have a 5 × 5 rotation matrix Â(α, β, γ) corresponding to rotations about the z-, x-and z-axes by the Euler angles α, β and γ respectively. The matrix elements a ij of  in terms of the matrix elements R ij of R are derived by writing the rotated real spherical harmonics, Y m in terms of the rotated Cartesian coordinates (x , y , z ) and then rewriting these in terms of a superposition of the original real spherical harmonics Y m . To demonstrate this, we provide one example of obtaining the transformation for m = -2. The rest of the rotated real spherical harmonics are derived in Appendix B.

Consider a generalized 3 × 3 Cartesian rotation matrix R which then gives us,

The matrix elements a ij are derived in Appendix B. Two solutions for matrix elements of the form a i3 are obtained. This is due to the fact that we have six quadratic terms (x 2 , y 2 , z 2 , xy, yz, zx) and only five coefficients a mm for a given m. However, we have shown that these two solutions are identical for rotations along each Cartesian axis and since every general rotation can be represented as rotations along the different axes, these two solutions are also identical for any general rotation.

The first step of our implementation is determination of the crystal-field geometry, since the final positions after rotation depend on the geometry (see Fig. 5). To distinguish between the octahedral and tetrahedral geometries, we first find the nearest neighbors and sort them by their distances. If the distance of the furthest away neighbor is over 35% of the closest, the crystal-field geometry is most likely tetrahedral. If this condition is not met, we determine the angles between the 'arms' of the octahedron. If these angles are within 35% of π/2, then the crystal-field is most likely octahedral. The reason we keep the relative angle and distance conditions to 35% is to include highly distorted octahedral configurations like those in GFO.

Once the crystal-field geometries are determined, we perform appropriate rotations about the z-, x-and z-axes by the Euler angles α, β and γ by setting the z-axis along along the direction where the angle formed between the central cation and its oxygen neighbors is the closest to π. These angles are obtained by using trigonometric relations between the final and initial positions. The final positions are 'known' since the neighbors of the cation have an arrangement as shown in Fig. 5. For an octahedral field, the six neighbors lie on the arms of the new coordinate system with the cation at the origin and their positions are along (±1, 0, 0), (0, ±1, 0) and (0, 0, ±1). In the case of a tetrahedral field, the four neighbors lie at the four vertexes of the two long diagonals in a cube with the cation at its center and the sides parallel to the axes of the new coordinate system. The positions are along

For a distorted octahedron, once the z-axis is set, as indicated above, the x and y axes are set by choosing the average angles between the cation and its two by two oxygen neighbors that approach best π/2.

For rotating the locally site projected wave functions, only the Euler angles to reach Thus, the rotation matrix for the new projector coefficients is the complex conjugate of the  rotation matrix. However, we are dealing with the real spherical harmonics, whose rotation matrices are completely real (see Appendix B). This implies that both the new projector coefficients and the real spherical harmonics transform using the same rotation matrix Â.

The implementation was done with the intention of studying the splitting of the e g and t 2g orbitals of a 3d cation in an octahedral or a tetrahedral environment. It is capable of handling distortions of up to 35 % in terms of the arm lengths or angles, as well as handling nonspin-polarized, spin-polarized and relativistic-spin-polarized electron wave functions. Our implementation is available as a Fortran90 module for VASP users. Let us assume that the real spherical harmonics for a given can be expressed in terms of a rotated real spherical harmonics set for the same . That is,

where we denote the rotation matrix as A = (a ij ), we have a 5×5 matrix for = 2. If we write out the full form of the spherical harmonics using Eq. A2 and using r ) 2 = r 2 , by multiplying out this matrix and expanding the primed coordinates in terms of unprimed (i.e. by using equation (B1)), we can compare the coefficients on the right and left sides of the equation to represent a ij in terms of R mn . Note that we have six polynomials (x 2 ,y 2 ,z 2 ,xy,yz,zx) but we wish to represent only 5 coefficients. This implies that we have a double solution for one coefficient. Indeed, if we follow through the calculation, we get two solutions for every a i3 .

Below is the list of all the elements of matrix A: