# Origin of magnetic frustration in BiMnO(No)

###### Abstract

BiMnO(NO) (BMNO) is a honeycomb bilayers anti-ferromagnet, not showing any ordering down to very low temperatures despite having a relatively large Curie-Weiss temperature. Using ab initio density functional theory, we extract an effective spin Hamiltonian for this compound. The proposed spin Hamiltonian consists of anti-ferrimagnetic Heisenberg terms with coupling constants ranging up to third intra-layer and fourth inter-layer neighbors. Performing Monte Carlo simulation, we obtain the temperature dependence of magnetic susceptibility and so the Curie-Weiss temperature and find the coupling constants which best matches with the experimental value. We discover that depending on the strength of the interlayer exchange couplings, two collinear spin configurations compete with each other in this system. Both states have in plane Néel character, however, at small interlayer coupling spin directions in the two layers are antiparallel (N state) and discontinuously transform to parallel (N state) by enlarging the interlayer couplings at a first order transition point. Classical Monte Carlo simulation and density matrix renormalization group calculations confirm that exchange couplings obtained for BMNO are in such a way that put this material at the phase boundary of a first order phase transition, where the trading between these two collinear spin states prevents it from setting in a magnetically ordered state.

###### pacs:

71.15.Mb, 75.50.Ee, 75.40.Mg.BiMnO(NO) (BMNO) is an experimental realization of frustrated honeycomb magnetic materials, synthesized by Smirnova et al. Smirnova et al. (2009). In this compound, the magnetic lattice can be effectively described by a weakly coupled honeycomb bilayers of Mn ions (Fig. 1). The temperature dependence of magnetic susceptibility of BMNO does not indicate any ordering down to K, in spite of the Curie-Weiss temperature K Smirnova et al. (2009); Onishi et al. (2012). The absence of long-range ordering in BMNO is also confirmed by specific heat measurements Smirnova et al. (2009); Onishi et al. (2012), neutron scattering Matsuda et al. (2010) and high-field electron spin relaxation (ESR) experiments Okubo et al. (2010). So far, the theoretical attempts to explain the magnetic properties of BMNO have been focusing on the frustration effect of second intra-layer coupling or the tendency toward dimerization by considering a large anti-ferromagnetic inter-layer nearest neighbor coupling Kandpal and van den Brink (2011); Ganesh et al. (2011a, b); Oitmaa and Singh (2012); Okubo et al. (2012); Zhang et al. (2014); Albarracín and Rosales (2016); Bishop and Li (2016). In an attempt to calculate the exchange interactions by ab initio method, it is found that the dominant exchange interactions are nearest in-plane coupling () and also an effective inter-plane coupling () which exceeds Kandpal and van den Brink (2011). However, in that work the experimental positions of the atoms in the structure are taken without geometry optimization and only magnetic configurations are considered for the calculation of exchange interactions. In this paper, we obtain a Heisenberg spin Hamiltonian for BMNO, using an ab initio LDA+U calculation. In our calculations, we consider a detailed analysis of nonidentical Mn atoms which were assumed to be identical in the previous calculation Kandpal and van den Brink (2011). We show that how this consideration can affect the exchange couplings in the spin Hamiltonian. We find that in contrary to the previous works, none of and are large enough to frustrate BMNO to reach an ordered state. Indeed, surprisingly the interlayer coupling constants are fined tuned in a way that make this system living at the edge boundary of two competing magnetic states.

Ab initio method. To derive magnetic exchange couplings, we employ Density Functional Theory (DFT) with Full-Potential Local-Orbital minimum-basis (FPLO), using FPLO code Koepernik and Eschrig (1999) (FPLO14.00-45). For charge analysis we employ Projector Augmented Wave (PAW) method, using Quantum-Espresso (QE) distribution Giannozzi et al. (2009). To account for exchange-correlation interaction we use PBE functional Perdew et al. (1996) from Generalized Gradient Approximation (GGA). To improve estimation of electron-electron Coulomb interactions, we also add Hubbard-like correction to DFT calculations, i.e., DFT+ Anisimov et al. (1991, 1993). To implement DFT+, FPLO uses Liechtenstein’s approach Liechtenstein et al. (1995); Eschrig et al. (2003). In Liechtenstein’s approach the two parameters, (on-site Coulomb repulsion) and (the on-site Hund exchange) needs to be set, which we use eV and and eV.

Spin Hamiltonian. The strategy of finding an effective spin Hamiltonian from ab initio calculations is to first compute the ground state energy for some given magnetic configurations. Then, mapping the energy difference of these configurations to an appropriate spin model gives us the coupling constants of the model. In this work, we use non-relativistic DFT, hence any magnetic anisotropy originating from the spin-orbit interaction is ignored in this approximation. Therefore, to leading order, we propose a spin Hamiltonian containing only bi-linear Heisenberg interactions, , where and are classical unit vectors representing the orientation of the magnetic moments at sites and , respectively, with exchange interactions between them. The primitive cell of BiMnO(NO) contains atoms. Because of the limitation in computational resources, we use the supercell containing atoms (Fig. 1). This lets us calculate ’s up to the third in-plane neighbor (, and ) and up to the fourth inter-plane neighbor coupling (, , and ) (see figure 2). BMNO is metallic in GGA, however, implementing spin-polarized calculation makes this compound insulating, independent of its magnetic configuration. Within co-linear spin polarized GGA, the ground state is a Néel state in which the nearest neighbor Mn magnetic moments (in and out of plane) are anti-parallel with respect to each other. This magnetic configuration is marked by N in Fig. 3. We calculated the total energy for more than independent magnetic configurations. Then employing the least square method, enables us to obtain the exchange couplings with the accuracy of meV. The top panel of figure 3 represents the energy landscape (per Mn atom) calculated for different magnetic configurations within the super-cell shown in Fig. 1. The detailed description of these configurations is given in Ref.sup (). As it is obvious from this figure, the two configurations N and N (bottom panel of Fig. 3) are very close in energy space. In configuration N, the magnetic ordering in each honeycomb layer is Néel type, but unlike N, the magnetic moment orientations of two layer are parallel. The coupling constants of the Heisenberg Hamiltonian, obtained by different values of on-site Coulomb interaction , are given in Table 1. We also checked that the exchange interactions between the adjacent honeycomb bilayers are negligible comparing to the ones inside the bilayers sup (). It is important to mention that to achieve equal couplings between equivalent Mn ions in the two layers, we need to geometrically optimize the atomic positions rather than just using the experimental atomic positions (for the details see supplementary information sup ()).

method | |||||||||
---|---|---|---|---|---|---|---|---|---|

FPLO | 1.5 | 10.7 | 0.9 | 1.2 | 3.0 | 1.1 | 0.5 | 0.9 | -244 |

2.0 | 9.0 | 0.8 | 1.0 | 2.6 | 0.9 | 0.5 | 0.8 | -203 | |

3.0 | 6.6 | 0.6 | 0.8 | 2.1 | 0.7 | 0.3 | 0.6 | -144 | |

4.0 | 5.1 | 0.5 | 0.6 | 1.7 | 0.6 | 0.3 | 0.5 | -111 |

The Mn spin state. The bond valence sum indicates the valence state BiMnO(NO3) for BMNO Smirnova et al. (2009). However, using the charge analyzing code Critic2 de-la Roza et al. (2014, 2009), within GGA/PAW, we find the valence state BiMnO(NO3) in N-configuration. This charge distribution will not change dramatically in the case of implementing DFT+ even with large U parameter. We also made sure that using an all-electron method, such as FPLO, this picture of charge distribution remains nearly unchanged. The local density analysis (Lowdin charges), also proposes the charge distribution BiMnO(NO3). These charge analyses show that the Mn-O bonds are ionic-covalent instead of being completely ionic. Indeed, the reason for such a fractional charge distribution in BMNO is the strong hybridization between Mn d-orbitals and the neighboring O p-orbitals. This also lowers the magnetic moment of Mn ions from to about (see Table I in supplementary information sup ()).

Monte Carlo Simulations. To gain insight into the finite temperature properties of the model Hamiltonian, we perform MC simulation. MC simulations are done on a supercell contains Mn atoms with periodic boundary conditions. We use single spin Metropolis updating, MC steps for thermalization and MC samplings for the measurement of physical quantity. To reduce the correlations, we skip MC sweeps between successive data collections. In Figure 4 the temperature dependence of inverse magnetic susceptibility for the sets of exchange couplings obtained by ab initio method using eV and also the experimental values are compared. This figure shows a very good agreement between the experimental values and the set of exchanges obtained by using eV. The linear fit at high temperatures crosses the -axis at a negative value which is the Curie-Weiss temperature . It can be seen that increases by increasing the value of onsite Coulomb repulsion . The is closest to what has been measured experimentally. To speculate about the ground state of the Hamiltonian, we calculated the inter-layer spin-spin correlation at a low temperature, up to fourth neighbor (Fig. 5). The correlations are calculated by averaging over MC samplings at . As it can be seen from this figure, for the couplings corresponding to eV (Table I), the spin-spin correlations between the two layers are very small. We observe that, increasing (decreasing) the value of interlayer coupling by a little amount pushes the system toward N (N) type ordering. Figure 5 shows the change of spin-spin correlation patterns, when varies by only percent, while keeping the rest of the couplings unchanged. The same results are obtained when only is changed while the other couplings are fixed (Fig.6) and also in the case that all the interlayer couplings are shifted up or downward (Fig.7). We also found that the small changes in the in-plane couplings do not induce spin ordering in the system.

Quantum effects. To make an inquiry about the quantum correlations at zero temperature, we use the density matrix renormalization group (DMRG) technique based on a matrix product state representation to evaluate the spin correlations functions Dolfi et al. (2014). In our calculations, we adopt (which is the spin of Mn) and the lattices with sites with . The spin-spin correlations normalized by (shown in the bottom panel of Fig. 5), confirms the transition from N to N states at , despite the weakening of the correlations as the effect of quantum fluctuations. This is while, the in-plane spins are in the Néel state, independent of the value of (Fig. 8).

Now, we proceed to determine the order of N-N transition for the set of exchange interactions obtained by using eV. Using Feynman-Hellmann theorem, the first derivative of the ground state energy with respect to a control parameter, say , is given by

(1) |

in which is the number of lattice points in a bilayer honeycomb of linear size and denotes the spin-spin correlation between the second interlayer neighbors. Figure 9 represents the DMRG results for the variation of this correlation versus , for the lattices of linear sizes . This figure shows a discontinuity in the second interlayer neighbor spin-spin correlation which becomes more pronounced by increasing the size of the system. Therefore equation 1 implies that the phase transition between these two ordered states is first order. Indeed, this result was expected because of different symmetries of N and N states.

Conclusion. In summary, we employed an ab initio LDA+U method to obtain the exchange coupling constants of a spin Hamiltonian for describing the magnetic properties of the honeycomb bilayer BMNO and figure out the reason that this compound does not show any ordering down to very low temperatures. Using eV, we found that a Hamiltonian containing only bilinear Heisenberg terms up to third in-plane and fourth out of plane neighbor, well matches the measured DC magnetic susceptibility for this material. Classical MC simulations and DMRG calculations on this spin Hamiltonian shows no sign of long-range ordering down to zero temperature. It is surprising that in BMNO the interlayer couplings are tuned in such a way to let this compound living at the phase boundary of the two collinear magnetic configurations N and N. Indeed, the interlayer coupling and encourage the N ordering, while and favor the N state. Therefore, the balance between these two sets of couplings adjusts BMNO to be at the N-N phase boundary. Hence, in the presence of any imbalance created as the effect of tension, compressive pressure, chemical doping, etc, the transition to N or N ordered states is expected. At this very special point, the spin-spin correlations in each layer are Néel type, however, there is almost a vanishing correlation between the two layers, making the dynamics of the two Néel states uncorrelated. The lack of correlations between the adjacent layer makes BMNO an effectively two-dimensional Heisenberg system for which there would be no finite temperature phase transition, according to the Mermin-Wagner theorem. It is also worthy to note that the presence of a strong enough spin-lattice interaction, could induce a spin-peierls lattice distortion and hence resolve the spin frustration. However, the reason that such a transition has not been observed experimentally, could be due to the small spin-orbit interaction in Mn atom which makes the temperature scale corresponding to the magneto-elastic interaction too small to cause an observable static lattice distortion in BMNO down to mK.

###### Acknowledgements.

We acknowledge Michel Gingras, Jeff Rau and Stefano de Gironcoli for the most useful discussions and comments. We also thank Phivos Mavropoulos for providing us with his MC code. DMRG results were checked by using the ALPS mps-optim application.## References

- Smirnova et al. (2009) O. Smirnova, M. Azuma, N. Kumada, Y. Kusano, M. Matsuda, Y. Shimakawa, T. Takei, Y. Yonesaki, and N. Kinomura, Journal of the American Chemical Society 131, 8313 (2009).
- Onishi et al. (2012) N. Onishi, K. Oka, M. Azuma, Y. Shimakawa, Y. Motome, T. Taniguchi, M. Hiraishi, M. Miyazaki, T. Masuda, A. Koda, et al., Physical Review B 85, 184412 (2012).
- Matsuda et al. (2010) M. Matsuda, M. Azuma, M. Tokunaga, Y. Shimakawa, and N. Kumada, Phys. Rev. Lett. 105, 187201 (2010).
- Okubo et al. (2010) S. Okubo, F. Elmasry, W. Zhang, M. Fujisawa, T. Sakurai, H. Ohta, M. Azuma, O. A. Sumirnova, and N. Kumada, Journal of Physics: Conference Series 200, 022042 (2010).
- Kandpal and van den Brink (2011) H. C. Kandpal and J. van den Brink, Physical Review B 83, 140412 (2011).
- Ganesh et al. (2011a) R. Ganesh, D. Sheng, Y.-J. Kim, and A. Paramekanti, Physical Review B 83, 144414 (2011a).
- Ganesh et al. (2011b) R. Ganesh, S. V. Isakov, and A. Paramekanti, Physical Review B 84, 214412 (2011b).
- Oitmaa and Singh (2012) J. Oitmaa and R. Singh, Physical Review B 85, 014428 (2012).
- Okubo et al. (2012) S. Okubo, T. Ueda, H. Ohta, W. Zhang, T. Sakurai, N. Onishi, M. Azuma, Y. Shimakawa, H. Nakano, and T. Sakai, Phys. Rev. B 86, 140401 (2012).
- Zhang et al. (2014) H. Zhang, M. Arlego, and C. Lamas, Physical Review B 89, 024403 (2014).
- Albarracín and Rosales (2016) F. G. Albarracín and H. Rosales, Physical Review B 93, 144413 (2016).
- Bishop and Li (2016) R. Bishop and P. Li, arXiv preprint arXiv:1611.03287 (2016).
- Koepernik and Eschrig (1999) K. Koepernik and H. Eschrig, Phys. Rev. B 59, 1743 (1999).
- Giannozzi et al. (2009) P. Giannozzi et al., Journal of Physics: Condensed Matter 21, 395502 (2009).
- Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Anisimov et al. (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
- Anisimov et al. (1993) V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993).
- Liechtenstein et al. (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
- Eschrig et al. (2003) H. Eschrig, K. Koepernik, and I. Chaplygin, Journal of Solid State Chemistry 176, 482 (2003), special issue on The Impact of Theoretical Methods on Solid-State Chemistry.
- (20) See Supplementary Material at [URL].
- de-la Roza et al. (2014) A. O. de-la Roza, E. R. Johnson, and V. LuaÃ±a, Computer Physics Communications 185, 1007 (2014).
- de-la Roza et al. (2009) A. O. de-la Roza, M. Blanco, A. M. PendÃ¡s, and V. LuaÃ±a, Computer Physics Communications 180, 157 (2009).
- Dolfi et al. (2014) M. Dolfi, B. Bauer, S. Keller, A. Kosenkov, T. Ewart, A. Kantian, T. Giamarchi, and M. Troyer, Computer Physics Communications 185, 3430 (2014).

## Supplemental Material

## Appendix A charge analysis

In this section, we explain why Mn-O bonds has ionic-covalent character. In fact, the ionic-covalent character of Mn-O bonds can be obsereved in the hybridization of Mn- and O- orbitals. The point group of MnO clusters in BiMnO(NO) is , where the -axis is perpendicular to honeycomb layer. Therefore, the Mn- orbitals are splitted as the effect of the crystal field into , and . The projected density of states plotted in figure 10 together with the orbital occupation calculation shown in Table. 2, indicate that the hybridization among orbitals of Mn and orbitals of O makes the crystal field states and to have fractional occupations.

Therefore the bond valence sums (BiMnO(NO3)) Smirnova et al. (2009) does not give rise to a true understanding of the charge distribution and also about magnetization of Mn atoms.

Figure 10 also shows a wide range of PDOS for all occupied -orbitals, which could be a reason that why the onsite electron-electron repulsion parameter is not very large in this system.

## Appendix B Exchange constants

In the experimental structure of BiMnO(NO) Smirnova et al. (2009) (with space group), within the experimental error, vertical positions of Mn atoms in the unit cell are almost equal (Table 3). In the DFT calculations, we found that Mn and Mn as well as Mn and Mn do not have exactly the same vertical positions, even if exact site geometry optimization is performed. However, as Table 3 shows, in DFT, the vertical positions of Mn and Mn and also Mn and Mn are very close. The site geometry optimization has been done within eV/Ang accuracy. The space group of geometry optimization structure and experimental structure is the same (P3, No 143). So the space group of Bi3Mn4O14 doesn’t change during geometry optimization.

Since Mn atoms are not completely identical, we have to use more Heisenberg constants. For example, instead of only one for first nearest neighbor interaction we need to consider two: one between and ( and the other between and (). For the second neighbor couplings we also have , , and . Similarly, there are variety of couplings for other inter and intra layer exchange interactions (see Table. 4).

To calculate the couplings of the Heisenberg Hamiltonian, we use magnetic configurations listed in Table. 6. Employing the least square method by considering these magnetic configurations, enables us to calculate the magnetic exchanges the within accuracy of 0.02 meV.

The ab initio results for the exchange couplings are given in Table. 4. These results are obtained after performing geometry optimization in a ferromagnetic configuration. The differences between the couplings of the same range are small (meV), hence their arithmetic mean are reported in the main paper.

Using the experimental structure, the difference between the couplings of the same range are significant (see Table 5). For example, =27.4 meV =10.8 meV within GGA/FPLO method. So assuming identical Mn atoms for experimental structure to derive are completely wrong.

In the main paper we assumed that the interaction between successive double layers could be ignored. To verify this assumption we calculate the energy difference of two magnetic configurations, a ferromagnetic configuration in which all the spins in all the bilayers are parallel and an anti-ferromagnetic configuration in which the spins in each bilayer are parallel, while they are antiparallel to the ones in neighboring bilayers. This energy difference gives an estimate of the strength of interaction between separated bilayers. Using FPLO method with eV, we obtain an energy difference (per Mn atom) of meV for these two configurations. This value is an order of magnitude less than the minimum energy difference (per Mn atom) of magnetic configurations inside each bilayer with respect to the ferromagnetic reference state, which is meV (Fig. 3 of main paper). Therefore it would be safe to neglect the bilayer-bilayer interaction.

method | |||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

FPLO | 1.5 | 10.8 | 10.7 | 1.0 | 0.8 | 1.0 | 0.8 | 1.2 | 1.2 | 2.9 | 3.0 | 1.1 | 1.1 | 0.5 | 0.6 | 0.9 | 0.9 |

2.0 | 9.0 | 9.0 | 0.9 | 0.7 | 0.9 | 0.7 | 1.0 | 1.0 | 2.6 | 2.6 | 1.0 | 0.9 | 0.4 | 0.5 | 0.8 | 0.8 | |

3.0 | 6.6 | 6.6 | 0.7 | 0.6 | 0.7 | 0.6 | 0.8 | 0.8 | 2.0 | 2.1 | 0.7 | 0.7 | 0.3 | 0.4 | 0.6 | 0.6 | |

4.0 | 5.1 | 5.1 | 0.5 | 0.4 | 0.5 | 0.5 | 0.6 | 0.7 | 1.6 | 1.7 | 0.6 | 0.6 | 0.3 | 0.3 | 0.5 | 0.5 |

method | |||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

FPLO | 0.0 | 27.4 | 10.8 | 0.9 | 0.7 | 2.5 | 2.0 | 1.2 | 2.2 | 4.6 | 5.9 | 1.8 | 1.3 | 0.2 | 0.6 | 1.0 | 0.7 |

2.0 | 18.4 | 6.5 | 0.6 | 0.5 | 1.6 | 1.4 | 0.8 | 1.3 | 2.7 | 3.6 | 1.2 | 1.0 | 0.3 | 0.6 | 0.9 | 0.7 | |

3.0 | 14.1 | 4.3 | 0.5 | 0.4 | 1.2 | 1.1 | 0.6 | 1.0 | 2.2 | 2.8 | 0.9 | 0.8 | 0.3 | 0.4 | 0.7 | 0.5 | |

4.0 | 11.1 | 2.8 | 0.3 | 0.3 | 1.0 | 0.8 | 0.5 | 0.8 | 1.7 | 2.3 | 0.7 | 0.6 | 0.2 | 0.3 | 0.5 | 0.4 |

## Appendix C Magnetic Configurations

In this section, we show the magnetic configurations used for the calculations of the Heisenberg coupling constants. To represent magnetic configurations, we assign a number on each Mn atoms in the supercell of BiMnO(NO) (figure 11), and then we specify the direction of Mn magnetic moments (up or down) by arrows shown in Table. 6.

Magnetic configuration | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn | Mn |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1 | ||||||||||||||||

2 | ||||||||||||||||

3 | ||||||||||||||||

4 | ||||||||||||||||

5 | ||||||||||||||||

6 | ||||||||||||||||

7 | ||||||||||||||||

8 | ||||||||||||||||

9 | ||||||||||||||||

10(N) | ||||||||||||||||

11 | ||||||||||||||||

12 | ||||||||||||||||

13 | ||||||||||||||||

14 | ||||||||||||||||

15 | ||||||||||||||||

16 | ||||||||||||||||

17 | ||||||||||||||||

18 | ||||||||||||||||

19 | ||||||||||||||||

20 | ||||||||||||||||

21 | ||||||||||||||||

22 | ||||||||||||||||

23 | ||||||||||||||||

24 | ||||||||||||||||

25 | ||||||||||||||||

26 | ||||||||||||||||

27 | ||||||||||||||||

28(N) | ||||||||||||||||

29 | ||||||||||||||||

30 | ||||||||||||||||

31 | ||||||||||||||||

32 | ||||||||||||||||

33 | ||||||||||||||||

34 | ||||||||||||||||

35 | ||||||||||||||||

36 | ||||||||||||||||

37 | ||||||||||||||||

38 | ||||||||||||||||

39 | ||||||||||||||||

40 | ||||||||||||||||

41 | ||||||||||||||||

42 | ||||||||||||||||

43 | ||||||||||||||||

44 | ||||||||||||||||

45 | ||||||||||||||||

46 | ||||||||||||||||

47(N | ||||||||||||||||

48 | ||||||||||||||||

49 | ||||||||||||||||

50 | ||||||||||||||||

51 | ||||||||||||||||

52 | ||||||||||||||||

53 | ||||||||||||||||

54 |

## Appendix D DMRG results analysis

We use density matrix renormalization group (DMRG) techniques based on a matrix product state (MPS) representation to evaluate the spin correlation functions. MPS is an variational ansatz that the variational parameters can be controlled by the matrix size, M, called the bond dimension. The ground state converges after sweeping a few times through the system.

We carry out 10 sweeps to converge the ground state within an error less than percent near the phase transition point. Moreover, we compare the correlation functions calculated using different bond dimension and find that the error of the results is less than percent. Away from the phase transition point, the errors are less than these values and reach to percent.

## References

- Smirnova et al. (2009) O. Smirnova, M. Azuma, N. Kumada, Y. Kusano, M. Matsuda, Y. Shimakawa, T. Takei, Y. Yonesaki, and N. Kinomura, Journal of the American Chemical Society 131, 8313 (2009).
- Onishi et al. (2012) N. Onishi, K. Oka, M. Azuma, Y. Shimakawa, Y. Motome, T. Taniguchi, M. Hiraishi, M. Miyazaki, T. Masuda, A. Koda, et al., Physical Review B 85, 184412 (2012).
- Matsuda et al. (2010) M. Matsuda, M. Azuma, M. Tokunaga, Y. Shimakawa, and N. Kumada, Phys. Rev. Lett. 105, 187201 (2010).
- Okubo et al. (2010) S. Okubo, F. Elmasry, W. Zhang, M. Fujisawa, T. Sakurai, H. Ohta, M. Azuma, O. A. Sumirnova, and N. Kumada, Journal of Physics: Conference Series 200, 022042 (2010).
- Kandpal and van den Brink (2011) H. C. Kandpal and J. van den Brink, Physical Review B 83, 140412 (2011).
- Ganesh et al. (2011a) R. Ganesh, D. Sheng, Y.-J. Kim, and A. Paramekanti, Physical Review B 83, 144414 (2011a).
- Ganesh et al. (2011b) R. Ganesh, S. V. Isakov, and A. Paramekanti, Physical Review B 84, 214412 (2011b).
- Oitmaa and Singh (2012) J. Oitmaa and R. Singh, Physical Review B 85, 014428 (2012).
- Okubo et al. (2012) S. Okubo, T. Ueda, H. Ohta, W. Zhang, T. Sakurai, N. Onishi, M. Azuma, Y. Shimakawa, H. Nakano, and T. Sakai, Phys. Rev. B 86, 140401 (2012).
- Zhang et al. (2014) H. Zhang, M. Arlego, and C. Lamas, Physical Review B 89, 024403 (2014).
- Albarracín and Rosales (2016) F. G. Albarracín and H. Rosales, Physical Review B 93, 144413 (2016).
- Bishop and Li (2016) R. Bishop and P. Li, arXiv preprint arXiv:1611.03287 (2016).
- Koepernik and Eschrig (1999) K. Koepernik and H. Eschrig, Phys. Rev. B 59, 1743 (1999).
- Giannozzi et al. (2009) P. Giannozzi et al., Journal of Physics: Condensed Matter 21, 395502 (2009).
- Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Anisimov et al. (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
- Anisimov et al. (1993) V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993).
- Liechtenstein et al. (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
- Eschrig et al. (2003) H. Eschrig, K. Koepernik, and I. Chaplygin, Journal of Solid State Chemistry 176, 482 (2003), special issue on The Impact of Theoretical Methods on Solid-State Chemistry.
- (43) See Supplementary Material at [URL].
- de-la Roza et al. (2014) A. O. de-la Roza, E. R. Johnson, and V. LuaÃ±a, Computer Physics Communications 185, 1007 (2014).
- de-la Roza et al. (2009) A. O. de-la Roza, M. Blanco, A. M. PendÃ¡s, and V. LuaÃ±a, Computer Physics Communications 180, 157 (2009).
- Dolfi et al. (2014) M. Dolfi, B. Bauer, S. Keller, A. Kosenkov, T. Ewart, A. Kantian, T. Giamarchi, and M. Troyer, Computer Physics Communications 185, 3430 (2014).