Orbital Ordering and Magnetic Interactions in BiMnO$_3$

This work is devoted to the analysis of orbital patterns and related to them interatomic magnetic interactions in centrosymmetric monoclinic structures of BiMnO$_3$, which have been recently determined experimentally. First, we set up an effective lattice fermion model for the manganese 3d bands and derive parameters of this model entirely from first-principles electronic structure calculations. Then, we solve this model in terms of the mean-field Hartree-Fock method and derive parameters of interatomic magnetic interactions between Mn ions. We argue that although nearest-neighbor interactions favors the ferromagnetism, they compete with longer-range antiferromagnetic (AFM) interactions, the existence of which is directly related with the peculiar geometry of the orbital ordering pattern realized in BiMnO$_3$ below 474 K. These AFM interactions favor an AFM phase, which breaks the inversion symmetry. The formation of the AFM phase is assisted by the orbital degrees of freedom, which tend to adjust the nearest-neighbor magnetic interactions in the direction, which further stabilizes this phase. We propose that the multiferroelectric behavior, observed in BiMnO$_3$, may be related with the emergence of the AFM phase under certain conditions.


Introduction
The multiferroic compounds have recently drawn an enormous attention due to promising practical applications as well as the fundamental interest [1,2,3,4]. Such systems, where magnetism coexists with the ferroelectricity, could be potentially used in the new devices aiming to transform the information in the form of the magnetization into the electric voltage and back. The primary goal of theorists is to unveil the microscopic mechanism leading to the coupling between magnetic and electric degrees of freedom.
The bismuth manganite (BiMnO 3 ), having highly distorted perovskite structure, has been regarded as one of the prominent multiferroic materials. Indeed, the ferromagnetism of BiMnO 3 is well established today. The Curie temperature (T C ) is about 100 K. The largest reported saturation magnetization is 3.92 µ B per one formula unit [5], 1 which is close to 4 µ B expected for the fully saturated ferromagnetic (FM) state. Nevertheless, the saturation magnetization decreases rapidly with the doping in Bi 1−x Sr x MnO 3 [6], that may indicate at the proximity of yet another and apparently antiferromagnetic (AFM) phase.
However, the situation around the ferroelectric properties of BiMnO 3 is more controversial. There are several facts, which do support the idea that BiMnO 3 is not only ferromagnetic, but also a ferroelectric material.
(i) The existence of ferroelectricity has been advocated by first-principles electronic structure calculations [7], and attributed to the chemical activity of the Bi(6s 2 ) lone pairs [8], in an analogy with other ferroelectric materials, such as PbTiO 3 .
(ii) According to early experimental data from electron and neutron powder diffraction, BiMnO 3 was considered to have noncentrosymmetric C2 space group in the entire monoclinic region [9,10], which is consistent with the ferroelectric behavior. Namely, BiMnO 3 undergoes two phase transitions at the temperatures of 474 and 770 K [9,10]. The first one at 474 K takes place without changing the monoclinic symmetry [11]. The phase transition at 770 K is monoclinic to orthorhombic [11] and was believed likely to be ferroelectric-paraelectric. Nevertheless, it is also worth to note that this point of view is rather controversial and according to [12] the onset of the ferroelectric behavior is expected only around 450 K, which in [12] was the point of isostructural (i.e., monoclinic to monoclinic) phase transition.
(iii) The ferroelectric hysteresis loop has been also reported in polycrystalline and thin film samples of BiMnO 3 [12], although the measured ferroelectric polarization was small (about 0.043 µC/cm 2 at 200 K). The first principle calculations performed for the experimental noncentrosymmetric structure result in much higher polarization (about 0.52 µC/cm 2 ) [13].
(iv) Kimura et al. [14] observed the changes of the dielectric constant induced by the magnetic ordering as well as by the external magnetic field near T C ∼ 100 K, and attributed them to the multiferroic behavior of BiMnO 3 .
(v) Sharan et al. [15] observed the electric-filed-induced permanent changes in the second harmonic response from the BiMnO 3 thin film, and argued that these changes are consistent with the possible presence of ferroelectricity.
However, there is also a growing evidence against the intrinsic ferroelectric behavior of BiMnO 3 .
After careful analysis, they concluded that both monoclinic phases observed in BiMnO 3 below 770 K have centrosymmetric space group C2/c. If so, BiMnO 3 should be an antiferroelectric, rather than the ferroelectric material. This funding was further confirmed in the neutron powder diffraction experiments by Montanari et al. [16] who also concluded that the crystal structure of BiMnO 3 is better described by the C2/c group in the wide range of temperatures (10 ≤ T ≤ 295 K) and magnetic fields (0 ≤ H ≤ 10 T). It is also important to note that there are many objective difficulties in the determination of the crystal structure of BiMnO 3 , which are mainly related with the strong effect of nonstoichiometry [17]. 2 (ii) For the related compound BiScO 3 , both neutron powder and electron diffraction analysis result in the centrosymmetric C2/c space group [19].
(iii) The structure optimization performed by using modern methods of electronic structure calculations revealed that the noncentrosymmetric C2 structure, which has been reported earlier [9,10], inevitably converges to the new total energy minimum corresponding to the C2/c structure with zero net polarization [20,21].
The goal of this work is to study of the orbital ordering and corresponding to it interatomic magnetic interactions in the centrosymmetric structure of BiMnO 3 . For these purposes we construct an effective lattice fermion model and derive parameters of this model from first-principles electronic structure calculations. After solution of this model we calculate the interatomic magnetic interactions. We argue that the peculiar orbital ordering realized below 474 K gives rise to FM interactions between nearestneighbor spins which always compete with longer-range AFM interactions. We propose that the ferroelectric behavior of BiMnO 3 can be related with the emergence of an AFM phase, which is stabilized by these longer-range interactions and breaks the inversion symmetry.
Thus, according to our point of view, the multiferroic behavior of BiMnO 3 is feasible and can be related with competition of two magnetic phases coexisting in a narrow energy range. The centrosymmetric FM ground state itself is antiferroelectric. Nevertheless, the ferroelectricity can be observed in the noncentrosymmetric AFM phase, which can apparently exist under certain conditions. Since the ferromagnetic (antiferroelectric) and antiferromagnetic (ferroelectric) phases can be stabilized by applying the magnetic and electric field, respectively, the magnetic moment can be switched off by the electric field and vice versa. This constitutes our idea of multiferroic behavior of BiMnO 3 . We rationalize several experimental facts on the basis of this picture.
The paper is organized as follows. In the next section we discuss details of the centrosymmetric crystal structure of BiMnO 3 . Section 3 briefly describes results of first principle electronic structure calculations in the local-density approximation (LDA). The construction of the model Hamiltonian is addressed in Section 4. The solution of the model Hamiltonian and physical meaning of interatomic magnetic interactions is discussed in Section 5. The results of calculations are discussed in Section 6. Finally, in Section 7 we will summarize our work and discuss how our results are related with the observed experimental behavior of BiMnO 3 .

Crystal Structure and Symmetry Considerations
BiMnO 3 has a highly distorted perovskite structure (figure 1). Figure 1. Fragment of crystal structure of BiMnO 3 , presented in the form of highly distorted perovskite lattice. The Bi atoms are indicated by the big yellow (light grey) spheres, the Mn atoms are indicated by the medium red (dark grey) spheres, and the oxygen atoms are indicated by the small hatched spheres. The primitive cell includes four Mn atoms, which are indicated by the numbers. The primitive translations are shown by arrows.
In our calculations we used the experimental crystal structure for T = 4 and 550 K obtained by Belik et al. The experimental structure parameters for T = 550 K can be found in [5], while the ones for T = 4 K are unpublished data [22]. 3 The primitive translations in the original monoclinic coordinate frame are give by 3 The structure parameters for T = 4 K are pretty close to the ones for reported in [5] for T = 300 K. If fact, all calculations have been performed using the experimental crystal structure both for T = 4 K and T = 300 K. Since both structures produce similar results, in the following we consider only the case of T = 4 K.
The space group C2/c has four symmetry operations: where in the notation {O|t}, O= E, I, m y , or C 2 y denotes the local symmetry operation, which is combined with the translation t= 0 or a 3 /2. Other notations are the following: E is the unity operation, I is the inversion, m y is the mirror reflection of the axis y, and C 2 y is the 180 • rotation around y.
The ferromagnetic configuration of BiMnO 3 has the full C2/c symmetry, which excludes any ferroelectricity. Possible antiferromagnetic configurations can be obtained by combining the symmetry operations (1) with the time-inversionT , which flips directions of the magnetic moments. Then, one can expect the following possibilities: (i) The AFM configuration ↑↑↓↓ (in these notations, four arrows correspond to the directions of the magnetic moments at the Mn-sites 1, 2, 3, and 4, respectively), which transforms according to the original space group C2/c. Similar to the FM case, this configuration exclude the ferroelectricity.
(ii) The AFM configuration ↑↓↓↑, which apart from the unity element {E|0}, has only one symmetry operation:T ⊗ {m y |a 3 /2}. This configuration does allow for the ferroelectricity, and the spontaneous polarization is expected to be perpendicular to the y-axis. Once the symmetry is broken by the AFM order, the atomic position will shift in order to minimize the total energy via magneto-elastic interactions. In this case, the crystal symmetry is expected to be P 2, which is compatible with the magnetic symmetry of the AFM ↑↓↓↑ phase. 4 Thus, the ferroelectric behavior in the ↑↓↓↑ phase is driven by the magnetic breaking of the inversion symmetry. A similar scenario of appearance of the ferroelectricity has been recently considered for other manganese oxides: HoMnO 3 [23,24] and TbMn 2 O 5 [25].
Other combinations ofT with the symmetry operations (1) will lead to unphysical solutions, where the local magnetic moments will vanish in one of the Mn-sublattices. Although such configurations are formally allowed by the symmetry, they clearly conflict with intraatomic Hund's first rule and are expected to have much higher energy. 5 We have also considered two ferrimagnetic configurations ↑↓↓↓ and ↓↓↑↓, which do not have any symmetry. In this case, the spontaneous polarization may have an arbitrary direction.

Electronic Structure in the Local-Density Approximation
First, we calculate the electronic structure corresponding to the low and high temperature structure of BiMnO 3 in the local density approximation (LDA) by using linear-muffin-tin-orbital (LMTO) method [26,27,28]. The atomic spheres radii and some details of the LMTO basis set used in the calculation are given in table 2. In order to fill the unit cell volume and reduce the overlap between atomic spheres, we additionally introduced 36 and 42 empty spheres for the 4 K and 550 K structure, respectively. The resulting total and partial densities of state are shown in figure 2.
The oxygen band, lying between -7 and -2 eV, is completely filled. The electronic structure near the Fermi level is mainly formed by the Mn(3d) states. Due to the hybridization, there is also a considerable weight of the Mn(3d) states in the oxygen band. The electronic structure near the Fermi level is further split into the Mn(e g ) and Mn(t 2g ) bands by pseudocubic crystal field operating in the MnO 6 octahedra, although in the highly distorted monoclinic structure there is no unique definition of the "t 2g " and "e g " orbitals since they are always mixed by the crystal distortion. The distortion is particularly strong in the low-temperature phase, leading even to an overlap between Mn(t 2g ) and Mn(e g ) bands. The Mn(e g ) band itself is split into the low-and highenergy subbands, lying at around 1 eV and 3 eV, respectively. There is a small gap between these subbands at around 1.7 eV. In total, the low-energy Mn(e g ) subband can accommodate one electron per one formula unit of BiMnO 3 . Therefore, according to the formal valence argument, in the fully polarized FM phase, both Mn(t 2g ) and low-energy Mn(e g ) bands are expected to be filled for the majority-spin channel, and the Fermi level is expected to fall in the pseudogap. The crystal distortion is somewhat released in the high-temperature structure, that leads to the opening of a gap between Mn(t 2g ) and Mn(e g ) bands and closing the gap between two Mn(e g ) subbands. The high-energy Mn(e g ) subbands always overlap with the Bi(6p) band spreading from about 2 to 5 eV. In this sense, any attempt to construct the model Hamiltonian for the isolated Mn(3d) bands will be an conjugated with some additional approximations for treating the Bi(6p) states and their hybridization with the Mn(3d) states.

Construction of the Model Hamiltonian
Our next goal is the construction of an effective multi-orbital Hubbard-type model for the Mn(3d) bands, located near the Fermi level, and derivation of the parameters of this model from the first-principles electronic structure calculations. The method has been proposed in [29]. Many details can be found in the recent review article [30].
The model itself is specified as follows: whereĉ † Rα (ĉ Rα ) creates (annihilates) an electron in the Wannier orbitalW α R centered at the Mn-site R, and α is a joint index, incorporating spin (s= ↑ or ↓) and orbital (m= xy, yz, z 2 , zx, or x 2 −y 2 ) degrees of freedom. The one-electron Hamiltoniant RR ′ = t α 1 α 2 RR ′ usually includes the following contributions: the site-diagonal part (R=R ′ ) describes the crystal-field splitting, whereas the off-diagonal part (R =R ′ ) stands for transfer integrals, describing the kinetic energy of electrons.
are the matrix elements of screened Coulomb interaction v scr (r, r ′ ), which are supposed to be diagonal with respect to the site indices {R}. The intersite matrix elements are typically small. Since we do not consider here the relativistic spin-orbit interaction, the matrix elements t α 1 α 2 RR ′ are diagonal with respect to the spin indices: i.e., t α 1 α 2 RR ′ = t m 1 m 2 RR ′ δ s 1 s 2 . The spin-dependence of the screened Coulomb interactions U R α 1 α 2 α 3 α 4 also has the regular form:

One-electron part
The one-electron part of the model Hamiltonian (2) can be constructed by using the formal downfolding method, which is applied to the Kohn-Sham equations within LDA. The method has been proposed in [29,31]. It is totally equivalent to calculation of the matrix elements of the Kohn-Sham Hamiltonian in the basis of Wannier functions constructed by using the projector-operator technique [32]. The advantage of the downfolding method is that by using it one can formally bypass the construction of the Wannier functions themselves and go directly to the calculation of the one-electron part of the model Hamiltonian. The comparison between original LDA bands and the ones obtained in the downfolding method is shown in figure 3. Generally, the agreement Figure 3. LDA energy bands corresponding to the low-temperature (left) and hightemperature (right) monoclinic structures of BiMnO 3 as obtained in the original electronic structure calculations using LMTO method (solid curves) and after the tightbinding (TB) parametrization using the downfolding method (dash-doted curves).
Notations of the high-symmetry points of the Brillouin zone are taken from [33].
is nearly perfect for the low-energy Mn(t 2g ) and the first four Mn(e g ) bands. In this region, the original electronic structure of the LMTO method is well reproduced after the downfolding. However, for the upper Mn(e g ) bands, which strongly overlap and interact with the Bi(6p) bands, it is virtually impossible to reproduce all details of the electronic structure in the minimal model consisting only of the Mn(3d) bands. 6 Therefore, in the upper-energy region, the electronic structure obtained in the downfolding method is only an approximation to the original LDA band structure. The model parameters for the one-electron part are obtained after the Fourier transformation of the downfolded Hamiltonian to the real space. In this case, the sitediagonal part oft RR ′ = t m 1 m 2 RR ′ describes the crystal-field splitting. The splitting of the e g levels in the low-temperature phase is particularly strong, being of the order of 1.5 eV (figure 4). As we will see below, the crystal-field effects will lead to a peculiar type of the orbital ordering, which will be mainly responsible for the magnetic properties of BiMnO 3 . This crystal-field splitting is mainly related with the change of the hybridization (or the covalent mixing) in different bonds of the distorted perovskite structure, which after the downfolding gives rise to the site-diagonal contributions in the model Hamiltonian. The nonsphericity of the Madelung potential, which plays a crucial role in the t 2g compounds [34,35], is considerably smaller then the effects of the covalent mixing in the e g bands and can be neglected. In the high-temperature phase, the crystal-field splitting shrinks  only in one of the sublattices, formed by the Mn atoms '1' and '2' in figure 1. In the second sublattice, formed by the atoms '3' and '4', the e g -level splitting remains large, being of the order 1 eV.
Because of complexity of the transfer integrals in the monoclinic structure, it is rather difficult to discuss the behavior of individual matrix elements of t m 1 m 2 RR ′ . Nevertheless, some useful information can be obtained from the analysis of averaged parameterst where d is the distance between Mn-sites R and R ′ . All transfer integrals are localized and practically restricted by the nearest neighbors at around 4Å (figure 5). The longerrange interactions are considerably smaller.

Screened Coulomb interactions
The matrix elements of screened Coulomb interactions in the Mn(3d) band,Û R , can be computed in two steps [29,30]. First, we perform the conventional constrained LDA calculations [36], and derive parameters of the on-site Coulomb interaction u = 10 eV and the intraatomic exchange interaction j = 1 eV. These parameters are certainly too large and only weekly depend on the crystal environment of Mn-atoms in the solid. They include several important channel of screening: for example, the screening of 3dinteractions by other electrons and the screening caused by relaxation of the atomic 3d wavefunctions are already included in the definition of u and j. However, this is not the whole screening and what the constrained LDA typically cannot do is to treat the socalled self-screening caused by the same 3d electrons, which participate in the formation of other bands due to the hybridization effects [30]. The major contribution comes from the O(2p) band, which has a large weight of the Mn(3d) states (figure 2). This channel . of screening can be efficiently taken into account in the random-phase approximation (RPA) by starting from interaction parameters obtained in the constrained LDA [29]: whereP is the static polarization matrix in RPA, calculated in the basis of atomic 3d orbitals, 7 andû is the 5×5×5×5 matrix of Coulomb interaction in the atomic limit. For each transition-metal site,û can be obtained from the parameters u and j by using a regular procedure, which is typically adopted in the LDA+U method [30]. The polarization matrixP is computed by using the LDA band structure. Nevertheless, in order to simulate the electronic structure close to the saturated ferromagnetic ground state, we have used different Fermi levels for the majority-and minority-spin channels in the process of calculation ofP . Namely, for the minority-spin channel, it was assumed that the Mn(3d) band is empty and the Fermi level has been placed right after the O(2p) band (i.e., around -1.5 eV in figure 3), while for the majority-spin channel, it was assumed that the Mn(3d) band accommodates all sixteen electrons (per four formula units). Meanwhile, we switch off all contributions to the polarization matrix related with the transitions between Mn(3d) bands in order to get rid of the unphysical metallic screening, which is present in RPA if one starts from the LDA band structure [37]. Then, for each Mn-site, we obtain the 5×5×5×5 matrixÛ R , which can generally depend on R and incorporate some effects of the local environment in solids. The symmetry of these matrices is also different from the spherical one. Nevertheless, just for the explanatory purposes, we fit each matrix in terms of three well known Table 3. Results of parametrization of screened Coulomb interactions in terms of the on-site Coulomb repulsion U , intraatomic exchange coupling J, and "nonsphericity" B for the low-temperature (T = 4 K) and high-temperature (T = 550 K) monoclinic phases of BiMnO 3 . All parameters are measured in eV. Positions of the Mn sites are explained in figure 1. parameters, which would fully specify all intraatomic interactions between 3d electrons in the spherical environment: the Coulomb repulsion U = F 0 , the intraatomic exchange coupling J = (F 2 + F 4 )/14, and the "nonsphericity" B = (9F 2 − 5F 4 )/441, where F 0 , F 2 , and F 4 are radial Slater's integrals. These parameters have the following meaning: U is responsible for the stability of certain atomic configuration with the given number of electrons, while J and B are responsible for the first and second Hund rule, respectively.
The results of such a fitting are shown in table 3. One can clear see that the Coulomb repulsion U is greatly reduced due to the self-screening effects, which are related with the admixture of the Mn(3d) states into the O(2p) band.

Hartree-Fock Approximation
In order to solve the model Hamiltonian (2) we employ the simplest mean-field Hartree-Fock approximation, where the trial many-electron wavefunction is searched in the form of a single Slater determinant |S{ϕ s k } , constructed from the one-electron orbitals {ϕ s k }. In this notation, k is a collective index combining the momentum k in the first Brillouin zone and the band number, and s is the spin of the particle. The one-electron orbitals are requested to minimize the total energy for a given number of particles N . The minimization is equivalent to the solution of Hartree-Fock equations for {ϕ s k }: wheret k ≡ t m 1 m 2 k is the one-electron part of the model Hamiltonian (2) in the reciprocal (similar equation for V ↓ Rm 1 m 2 is obtained by interchanging ↑ and ↓). Equations (4) are solved self-consistently together with the equation for the density matrixn s ≡ n s Rm 1 m 2 in the basis of Wannier functions. After self-consistency, the total energy (3) can be computed as

Magnetic Interactions
By knowing {ε k } and {ϕ k }, one can construct the one-electron (retarded) Green which can be used in many applications. For example, the interatomic magnetic interactions corresponding to infinitesimal rotations of the spin magnetic moments near the equilibrium can be easily computed as [38,39]: where ∆V R =V ↑ R −V ↓ R is the magnetic (spin) part of the Hartree-Fock potential, Tr L denotes the trace over the orbital indices, and ε F is the Fermi energy. According to (6), J RR ′ >0 (<0) means that for the given magnetic state, the spin alignment in the bond RR ′ corresponds to the local minimum (maximum) of the total energy. However, in the following we will use the universal notations, according to which J RR ′ >0 and <0 will stand the ferromagnetic and antiferromagnetic coupling, respectively. These notations correspond to the mapping of the total energy change of the Hartree-Fock method, associated with the small rotations of the magnetic moments, onto the Heisenberg model [38]: where e R is the direction of the spin magnetic moment at the site R. Generally, the parameters {J RR ′ } are not universal and depend on the magnetic state in which they are calculated (for example, through the change of the orbital ordering [35] or the electronic structure [40] in each magnetic state).
If we are dealing with the collinear magnetic structure, where all spins are parallel to the z-axis, i.e. e R = (0,0,1) or (0,0,−1), one can consider a small rotation of the magnetic moment at one of the site, e R = (cos θ R sin φ R , sin θ R sin φ R , cos θ R ), and calculate the second derivative of E Heis with respect to θ R : In this expression, s RR ′ = 1 and −1 stands correspondingly for the FM and AFM alignment in the bond RR ′ . J 0 R characterizes the stability of the magnetic system with respect to the rotation of the single spin. It can be also related with the spin stiffness and the magnetic transition temperature in the mean-field approximation [38]. If J 0 R > 0, the spin system is stable while if J 0 R < 0 it is unstable.

Decomposition into "Double Exchange" and "Superexchange"
Many properties of perovskite manganese oxides are related with the simple fact that the exchange splitting ∆V R is large, and for many applications can be treated as the largest parameter in the problem [40,41]. This is because Mn 3+ ions have four unpaired 3d electrons, which interact through the Hund's rule coupling J. Loosely speaking, the exchange splitting between the majority and minority spin states is controlled by the parameter U+3J, which is about 4.9 eV (table 3), whereas the orbital polarization (or the splitting of occupied states with one particular projection of spin) is controlled by U−J, being "only" about 1.4 eV. Therefore, as the first approximation, one can neglect the orbital dependence of ∆V R end replace it by some constant exchange splitting ∆ ex : i.e., A typical example of the exchange splitting in the low-temperature monoclinic phase is shown in figure 6: the averaged exchange splitting ∆ ex is about 4.7 eV, whereas the deviations from ∆ ex for the particular orbitals do not exceed 1.5 eV. Of course, (8) is a crude approximation. Nevertheless, as will see below, it appears to be very useful for the analysis of interatomic magnetic interactions. It also reproduces the main trends of the behavior of these interactions at least on the semi-quantitative level. Since ∆ ex is large, all minority-spin states are empty ( figure 7). Therefore, all poles ofĜ ↓ R ′ R are located in the unoccupied part of the spectrum and below ε F one can use the 1/∆ ex expansion forĜ ↓ R ′ R [40,41]. Then, the first two terms in the expansion of J RR ′ will have the following from: and where we have used the notationsĥ R ′ R =t R ′ R +V ↑ R δ R ′ R and (ĥ) 2 R ′ R = R ′′ĥR ′ R ′′ĥ R ′′ R . J D RR ′ is proportional to {t R ′ R } and does not depend on ∆ ex . In an analogy with [40,41], we will called it "the double exchange interaction", although, strictly speaking, it is not a regular double exchange sinceĜ ↑ RR ′ also includesV ↑ R , which takes into accounts the effects of the orbital polarization of the electronic origin. J S RR ′ incorporates the effects of the second order with respect to {t R ′ R } and is inversely proportional to ∆ ex . Therefore, in the following it will be called "the superexchange interaction". We will also consider two approximation for J D RR ′ and J S RR ′ . In the first one,ĥ R ′ R will be regular Hartree-Fock Hamiltonian for the majority-spin states andĜ ↑ RR ′ is the Green function corresponding to this Hamiltonian: In the second one, in order to be consistent with the approximate expression (8) for the exchange splitting, we will neglect all effects of the orbital polarization of the electronic origin also in the definition ofĥ R ′ R andĜ ↑ RR ′ . Therefore, apart from the constant shift, h R ′ R is replaced byt R ′ R . Then,Ĝ ↑ RR ′ is the regular LDA Green function, which is obtained from (11) after replacingĥ byt. In this approximation, it becomes more clear why we continue to use the term "double exchange", even though our system can be insulating, like the low-temperature phase of BiMnO 3 , where already in LDA there is a gap between Mn(e g ) bands (figure 3). It is true that the existence of this gap, ∆, is related with some kind of the orbital polarization. In this sense it is still reasonable to consider the superexchange processes by treating all transfer integrals as a perturbation. This would correspond to the superexchange interactions of the form t 2 eff /∆, where t 2 eff is the square of an effective transfer integral between Mn sites, which is related with {t R ′ R }. However, this orbital polarization comes from the large crystal-field splitting, which is just another effect of the covalent mixing, and, therefore has the same origin as t eff . Therefore, ∆ should be proportional to t eff , and the "superexchange" t 2 eff /∆ becomes also proportional to t eff . From this point of view, it is still reasonable to call this interaction as the "double exchange".

Results and Discussion
The orbital ordering in the low-temperature monoclinic phase of BiMnO 3 is shown in figure 8. We have tried three different methods in order to derive the distribution of the electron density around Mn-sites.
(i) In the first method, we simply calculate the site-diagonal elements of the density matrix in the original LMTO basis by integrating over four lowest e g bands (spreading around 1 eV in figure 3), and plot the electron density corresponding to this density matrix.
(ii) In the second method, we plot the densities of the lowest e g -orbitals obtained from the diagonalization of the site-diagonal part of the one-electron Hamiltonian, which was derived from the downfolding method. In the other words, these are just the crystal-field orbitals, corresponding to the fourth atomic level in figure 4.
(iii) In the third method, we plot the electron density for the occupied ↑-spin e g band, obtained in the Hartree-Fock calculations for the ferromagnetic state (figure 7).
All three methods provide a very consistent picture for the general details of the orbital ordering in the low-temperature phase of BiMnO 3 , which is also consistent with results of full-potential calculations by Shishidou [20]. The orbital ordering in the high-temperature structure, derived from the crystalfield orbitals, is shown in figure 9. This orbital ordering is entirely related with the crystal distortion, which splits the atomic e g levels. Even for the sites '1' and '2' with the least distorted environment (table 4), this splitting is of the order of 0.2 eV. This corresponds to the temperature of about 2100 K, which largely exceed the temperature of monoclinic-to-orthorhombic transition (about 770 K [14]). Thus, it is reasonable to expect the orbital ordering to take place in both monoclinic phases, below and above 474 K. However, as it is clearly seen from the comparison of figures 8 and 9, the character of the orbital ordering will change at the point of phase transition. This conclusion is qualitatively consistent with results of resonant x-ray scattering on BiMnO 3 [42].
Results of Hartree-Fock calculations of the total energies for the ferromagnetic and several antiferromagnetic configurations are shown in table 4. Our main observation is that for the low-temperature monoclinic structure the ferromagnetic phase appears to be nearly degenerate with the antiferromagnetic ↑↓↓↑ phase, which can be obtained by flipping the directions of the magnetic moments at the Mn-sites '2' and '3'. Perhaps, the tendencies towards the antiferromagnetism are somewhat overestimated in our model, and there are several reasons for it: (i) The calculations of the Coulomb interaction U are always conjugated with certain approximations [29,30], and some of these parameters may be underestimated. As we will see below, larger values of the parameter U would indeed help in stabilizing  (ii) Our model (2) does not explicitly include the oxygen states. This appears to be a good approximation for titanium and vanadium perovskite oxides [35], where the transition-metal and oxygen bands are well separated. However, manganese compounds are much closer to the charge-transfer regime because of the proximity of Mn(3d) and O(2p) bands, and much stronger hybridization, which takes place between these groups of states. Moreover, as it was pointed out in section 4.1, there is also an overlap between Mn(3d) and Bi(6p) bands. Therefore, the Hubbard model (2), where the form of the Coulomb interactions is borrowed from the atomic limit for the Mn(3d) states is an approximation, which may ignore some contributions to the relative stability of different magnetic configurations. For example, it is known that the magnetic polarization of the oxygen states will additionally stabilize the FM phase [43]. These effects are not included into the model (2).
Nevertheless, as we will see below, the competition between ferromagnetic and antiferromagnetic ↑↓↓↑ phases itself is a genuine effect, which is directly related with the form of the orbital ordering in the low-temperature monoclinic structure. The distance-dependence of interatomic magnetic interactions calculated in the lowtemperature monoclinic structure is shown in figure 10. These calculations have been performed using the formula (6) for infinitesimal rotations of magnetic moments near the ferromagnetic state. We note the following. There are two types of relatively strong ferromagnetic interactions between nearest neighbors, which operate in the bonds 1-3 and 1-3 ′′ (see figure 11 for notations). 8 The character of these interactions is directly related with the "antiferromagnetic" orbital ordering in the bonds 1-3 and 1-3 ′′ , and can be anticipated already from distribution of the Mn-O bondlengths in the lowtemperature monoclinic phase [5]. 9 The interaction in the third bond 1-3 ′ , formed by the Figure 11. Fragment of the orbital ordering pattern in the low-temperature monoclinic phase of BiMnO 3 . Different magnetic sublattices, which are formed in the antiferromagnetic ↑↓↓↑ structure are shown by different colors.
nearest neighbors, is relatively week. This is again consistent with the geometry of the orbital ordering, corresponding to the minimal overlap between occupied and unoccupied e g orbitals. Similar situation occurs in the bond 1-4 ′ . The most striking result of the present calculations is the existence of relatively strong long-range antiferromagnetic interaction in the bond 1-2 ( figure 11). Nevertheless, this result is also anticipated from the geometry of the orbital ordering. Note that the occupied e g orbitals at the sites 1 and 2 are directed towards each other. Although the direct transfer integrals between these two sites are relatively small (see figure 5), it is still reasonable to expect the existence of AFM interactions, which are mediated by unoccupied e g states of the site 3 ′′ . Such a situation is somewhat similar to superexchange interactions, which take place via oxygen states in the charge-transfer insulators [45,46], and the mechanism itself is sometimes called as the "super-superexchange". Thus, in the low-temperature monoclinic phase of BiMnO 3 we are always dealing with a competition of nearest-neighbor FM and longer-range AFM interactions. In fact, there are several factors, which can make these interactions comparable with each other. It is true that, generally, the nearest-neighbor interactions are expected to be much stronger, because all transfer integrals are basically restricted by the nearest neighbors ( figure 5). However, for the nearest-neighbor interactions we are also dealing with the strong cancelation of FM "double exchange" and AFM superexchange contributions (table 5). For example, this cancellation is nearly perfect for the "weak bonds" 1-3 ′ and therefore, favors the ferromagnetic interactions between these sites [44]. 1-4 ′ . This is a general rule for perovskite manganese oxides, which explains a strong reduction of nearest-neighbor magnetic interactions, so that they can easily become comparable with some longer-range interactions [39,40,41]. On the other hands, for the longer-range AFM interaction in the bond 1-2, there is no such cancellation. The longrange interactions are expected to vanish for the undoped (parent) manganites, provided that they would have a undistorted cubic structure. This effect is entirely related with the symmetric filling (or half-filling) of the majority-spin e g band [41]. Nevertheless, many parent manganites (like BiMnO 3 ) have a strongly distorted crystal structure. This distortion gives rise to the orbital ordering, which leads to certain asymmetry of filling of the majority-spin e g band, and this asymmetry is finally manifested in the appearance of longer-range interactions.
Since nearest-neighbor interactions favor the ferromagnetism, while the longerrange interactions favor the formation of the antiferromagnetic ↑↓↓↑ structure, one can generally expect the competition between these two phases, as it is clearly seen from results of total energy calculations shown in table 4. Nevertheless, there is another factor, which will additionally stabilizes the AFM ↑↓↓↑ phase and explains why it is so close in energy to the FM phase.
Note that the orbital degrees of freedom can additionally change their form in order to modify the interatomic magnetic interactions in the direction which will further stabilize the given magnetic structure [44]. Since the form of the orbital ordering is efficiently constrained by the large crystal-field splitting, it does not strongly depend on the magnetic state, and visually one can observe only tiny changes in the distribution of occupied e g electron density ( figure 12). Nevertheless, it is interesting to see that these tiny changes may have a profound effect on the behavior of interatomic magnetic interactions. Indeed, in the AFM ↑↓↓↑ structure, the chain 1-3 ′′ -2 contains one AFM bond (1-3 ′′ ) and one FM bond (2-3 ′′ , see figure 11). Although in the AFM structure, both interactions remain ferromagnetic, there is a clear polarity of interactions and the magnetic coupling in the AFM bond 1-3 ′′ is considerably weaker than the one in the ferromagnetic bond 2-3 ′′ (table 5). Even more dramatic change occurs in the plane of the distorted perovskite structure. In the AFM structure, even visually one can see some anisotropy in distribution of the occupied e g electron density at the site 1 ( figure  12), which appears to be more contracted in the direction of the FM bond 1-4. On the other hand, this distribution is nearly isotropic in the FM phase. This means that in the direction 1-4 of the AFM phase, the weight of the e g orbitals is additionally moved into the unoccupied part of the spectrum. This opens some additional pathes for the virtual hoppings into the unoccupied part of the spectrum, which will additionally stabilize the FM coupling. Indeed, the exchange coupling in the FM bond 1-4 is 7.8 meV, while the one in the AFM bond 1-3 is strongly reduced till 1.7 meV. Moreover, the magnetic interactions in the weak bonds 1-3 ′ and 1-4 ′ are also adjusted by the orbital-ordering effects. For example, the AFM coupling in the bond 1-3 ′ is enhanced, while the coupling in the FM bond 1-4 ′ becomes ferromagnetic.
Because of these orbital ordering effects both magnetic configurations appear to be locally stable. Indeed, the parameters J 0 R , calculated in the FM and AFM ↑↓↓↑ phases are (7.4, 12.1) and (17.6, 8.9) meV, respectively, where the first number in parenthesis corresponds to the Mn-site 1, while the second number corresponds to the Mn-site 3. All parameters are positive, meaning that both configurations are stable, at least with respect to independent rotations of magnetic moments at the sites 1 and 3.
In addition to electronic degrees of freedom, the AFM ↑↓↓↑ phase can be additionally stabilized by structural effects associated with polar atomic displacements in the direction which further minimizes the total energy of the system via magnetoelastic interactions. For example, magnetic couplings in the bonds 1-3 ′′ and 2-3 ′′ can be further adjusted by displacement of the Mn-atom 3 ′′ ( figure 11). Although we do not consider such a mechanism here (that would require detailed structural optimization using modern full-potential electronic structure methods), it would be certainly an interesting step to do in the future. This mechanism was considered for example in [25,24] for other multiferroic compounds.
In the high-temperature monoclinic structure, the FM phase is clearly the most stable one (table 4). This is closely related with the fact that the FM phase is metallic (figure 7) and the double exchange interactions clearly dominate (table 5). However, this is true only for the low-temperature regime. At elevated temperatures, the structure of interatomic magnetic interactions will be largely modified by the magnetic disorder, which may also destroy metallic character of the electronic structure. 10 Thus, it is rather meaningless to discuss the properties of the high-temperature phase in the lowtemperature limit. The magnetic disorder in the high-temperature phase is certainly one of the interesting problems. However, it is beyond the scopes of the present work. The possible tools to address this problem is the coherent potential approximation (CPA) [47,48], or, more generally, the dynamical mean-field theory (DMFT) [49]. Nevertheless, we would like to emphasize that since the type of the orbital ordering changes above 474 K (figure 9), the long-range interactions between atoms 1 and 2 are expected to be weak. Therefore, the long-range AFM correlations, which according to our point of view are indispensable for the inversion symmetry breaking and appearance of the ferroelectricity, should not play any important role in the high-temperature monoclinic phase. Thus, according to our scenario, the temperature of the isostructural phase transition (474 K) should be also regarded as a upper bound for the possible onset of ferroelectricity.

Implications to the Properties of BiMnO 3
Thus, we would like to propose that the multiferroic behavior of BiMnO 3 should be closely related with a competition between two magnetic phases. One is the centrosymmetric (and antiferroelectric) FM phase. Another one is the AFM ↑↓↓↑ phase, which breaks the inversion symmetry and allows for the spontaneous electric polarization in the direction perpendicular to the y-axis. The existence of both phases is closely related with the peculiar orbital ordering, which takes place below 474 K.
This means that despite the fact that BiMnO 3 is crystallized in the centrosymmetric C2/c structure, there is still a room for the multiferroic behavior, if we could engineer the samples where these two phases coexist in a narrow energy range accessible for the physical changes of electric and magnetic fields as well as the temperature T . Then, one could readily expect the "switching phenomena". For example, by applying the magnetic field H one could stabilize the FM antiferroelectric phase and switch off the net electric polarization. Conversely, one could apply the electric field E and stabilize the AFM ferroelectric phase with zero net magnetization.
It is true the low-temperature ferromagnetism in the pure bulk samples of BiMnO 3 is well established today [5]. Therefore, the symmetry is expected to be C2/c and these samples are not extremely promising from the viewpoint of multiferroic applications. However, it is also known that the BiMnO 3 is extremely difficult to synthesize, especially in the single crystalline form. Therefore, it is reasonable to expect that the BiMnO 3 samples will always have some defects, and we would like to speculate that these defects may play a positive role in stabilizing some fractions of antiferromagnetic and noncentrosymmetric phase. This seems to be reasonable, because the low-temperature magnetization observed in the BiMnO 3 samples, which do exhibit the ferroelectric behavior, was only 2.6µ B and reached 3.1µ B per one formula unit in high magnetic fields [12]. These values are considerably lower that 4µ B expected for the single saturated FM phase, meaning that the samples were not sufficiently pure and might contain a fraction of the AFM phase. This strongly reminds the idea of clustering or macroscopic phase separation, which has been intensively discussed in other perovskite manganese oxides, in the context of colossal magnetoresistance phenomena [50,51].
In this respect it is important to note that even in the high-quality samples, whose low-temperature saturation magnetization was close to 4µ B [5], the authors of [17] and [52], by means of atomic pair distribution function analysis on neutron powder diffraction data and selected-area electron diffraction technique, respectively, have observed the existence of short-range ordered structures (or domains) with the broken inversion symmetry. The symmetry of these structures was either P 2 or P 2 1 , which is consistent with the crystal symmetry of the AFM ↑↓↓↑ phase (P 2), that we propose. Moreover, these experimental data strongly suggests that the breaking of the inversion symmetry is mainly caused by the Mn atoms, that is again consistent with the idea of the magnetic origin of this effect. Unfortunately, the authors of [17,52] focused only on the structural properties of BiMnO 3 , and did not provide any information on how these structural properties can be related with the magnetic behavior of BiMnO 3 . From our point of view, such measurements would be very useful. For example, if the origin of the noncentrosymmetric domains was indeed magnetic, it is reasonable to expect the size and relative wight of these domains to decrease in the magnetic field.
Another possibility to control the properties of BiMnO 3 is to use the thin films. In this case, due to the lattice mismatch of the bulk BiMnO 3 and the substrate, the latter causes an additional strain and may strongly affect the magnetic behavior of BiMnO 3 . This effect is well known for other (colossal magnetoresistive) manganese oxides, and can be used as an efficient tool for controlling the electronic and magnetic properties of these systems near the point of phase transition between FM and AFM states [53,54]. Moreover, the magnetic structure at the surface of BiMnO 3 may be also different from the one in the bulk. Indeed, the saturation magnetization of the BiMnO 3 thin films grown on the (100) SrTiO 3 substrate was only 2.8µ B [55], which is considerably smaller than the bulk value. The magnetic moment increases almost linearly with the increases of the film thickness and reaches the nearly saturated value at around 500Å. The use of the (110) SrTiO 3 substrate yields even smaller saturation moment (about 1.8µ B [55]). Thus, all these data suggest that the magnetic ground state realized in the thin films of BiMnO 3 is not a pure FM one, and may contain some elements of the AFM structure. Of course, on the basis of this highly limited experimental information about the saturation magnetization it is impossible to make a ultimate conclusion whether this AFM structure is indeed the ↑↓↓↑ one, which we propose. Nevertheless, one can speculate that the ferroelectric behavior, which is apparently seen in the BiMnO 3 thin films, can be again related with the deviation of the magnetic structure from the pure FM one [12,15].
Finally, although the ground state of BiMnO 3 is ferromagnetic, some fraction of the AFM ↑↓↓↑ phase can emerge at elevated temperatures. If at zero temperature both FM and ↑↓↓↑ AFM phases correspond to the local minima of the total energy and are connected to each other by the first-order transition, the rising of the temperature can naturally lead to a coexistence of these two states. By neglecting the interaction between two phases, this effects is simply related with the configuration mixing entropy. Moreover, the appearance of the two-phase state is considerably facilitated in the presence of defects [50,51], as it is well known for other (colossal magnetoresistive) manganites [56]. This behavior implies the existence of the temperature hysteresis loop in the magnetization curve near T C . A small temperature hysteresis has been indeed reported in [5], which may be related with small anomalies of dielectric constant observed in [14]. Since the existence of the long-range AFM interaction is closely related with the peculiar orbital ordering persisting up to 474 K, no ferroelectricity can be generally expected above this temperature.
In summary, we believe that further exploration of multiferroic behavior of BiMnO 3 should be focused on the revealing of the AFM ↑↓↓↑ phase. Possible multiferroic applications of BiMnO 3 will strongly depend on whether one can find the conditions of coexistence of the FM and AFM ↑↓↓↑ phases in a narrow energy range accessible for the switching by electric and magnetic fields. It seems that the available experimental data do not rule out this idea, although at the present stage there is no direct support to it either. We hope that our work will stimulate theoretical and experimental activity in this direction.