Theoretical study on the role of dynamics on the unusual magnetic properties in MnBi

We study the electronic structure and lattice dynamics in the ferromagnet MnBi using first-principles calculations and a tight-binding model. The band structure around the Fermi level is dominated by Bi-p states which are the primary contributors to the magnetic anisotropy energy in the low temperature structure. A tight-binding model consisting of Mn-d and Bi-p states is developed and the parameters are determined from first-principles calculations. Phonon dispersions and elastic moduli exhibit several interesting features. The results imply that the magnetic interaction with the crystal lattice in MnBi is considerably more complex than previously thought and in particular that there is a rich interplay between phonons and magnetism involving both magnetoelastic and magnetostrictive coupling.

I n recent years, significant attention has been paid to realize permanent magnets without rare earth elements. Due to its large uniaxial magnetic anisotropy energy (MAE) and Curie temperature above room temperature, MnBi has emerged as a favorable candidate among transition metal systems [1][2][3] . In addition, it has good magneto-optical properties including a large Kerr rotation 4 with potential applications in erasable magnetooptical memory devices 5 .
In the low-temperature phase, MnBi is a ferromagnetic compound with the NiAs structure. Above 628 K, a first-order paramagnetic transition accompanied by a phase decomposition to Mn 1.08 Bi takes place 6 . The magnetic phase of MnBi is unusual in that the anisotropy constant (K u ) increases with temperature unlike most magnets 2 . Also, K u is negative below ,90 K, with all the magnetic moments aligned in the in-plane direction. A spin-reorientation transition takes place at T SR < 90 K, when the magnetization changes direction towards c axis. Thus, K u starts at 20.2 MJ/m 3 close to 0 K and reaches 2 MJ/m 3 at room temperature. The lattice parameters also exhibit a small kink at the spin-reorientation temperature indicative of a first order transition 7 . Mechanical grinding of the crystalline samples were found to enhance room temperature coercivity, at the same time suppressing the spin-reorientation temperature 8 which was attributed to the competition between second and higher order crystal field terms. MnBi also exhibits a large Kerr effect, with a double peak spectrum similar to other ferromagnetic transition metals 4,9 . Interband transitions between Mn 3d and Bi 6p states were suggested to be responsible for the effect.
Several theoretical studies have been carried out to understand the magnetic properties of MnBi 3,10-13 . Electronic structure calculations correctly reproduce the low temperature magnetic structure, which is ferromagnetic and metallic 14 . Calculated K u is negative, in agreement with experiments, but it's magnitude is overestimated. Sakuma et al found that a slight decrease in the valence electron number can change the sign of K u , so that partially substituting Bi with Sn (MnBi 12x Sn x ) can lead to the positive anisotropy observed at high temperatures 11 . The volume derivative of K u is also negative in density functional theory (DFT) based calculations, but including Coulomb correlations through DFT 1 U makes it positive 12,13 . Thus, thermal expansion of the lattice was suggested as a potential reason for the spin-reorientation transition by Zarkevich et al, who found MAE to change sign for larger values of the lattice parameter a 12 . More recently, Antropov et al 13 studied the spinreorientation in MnBi with the help of extensive LSDA 1 U calculations, in which they treated U as an adjustable parameter and found good agreement with experimental observations for a larger value of U. Since, MnBi is not expected to be a strongly correlated system, it is surprising that anisotropy energies depend strongly on U.
In MnBi, substantial coupling of magnetism with the lattice is implied, which may lead to interesting lattice dynamical behavior in relation to magnetism. Thermal vibrations of atoms can alter the crystal field of the ions which can in turn affect the single ion anisotropy. In addition, spin-fluctuations lead to temperature dependence of the molecular field, also affecting the anisotropy. These effects have been suggested as responsible for spin-reorientation in rare earth magnets such as Nd 2 Fe 14 B 15 . Hence, in this paper, we explore the role of lattice dynamics on the magnetic behavior of MnBi with the help of first principles calculations. We also present a simple tight-binding model for the low temperature structure involving 3d and 6p electrons, which can reproduce the magnetic ground state. We find that calculated phonon spectra and sound speeds in MnBi are high anisotropic. Our calculations suggest that thermally induced lattice vibrations contribute significantly to the spin-reorientation transition in addition to thermal expansion.

Results
Effect of Hubbard U. In the low temperature phase, MnBi crystallizes in the hexagonal structure with space group P 63 =mmc{D 4 6h , with the unit cell as shown in Fig. 1. The experimental structure from Ref. 6 is used in our calculations, with lattice constants a 5 b 5 4.285 Å and c 5 6.113 Å . Along the c direction, Mn and Bi are arranged in alternating layers, each bonded to six nearest neighbors. The Bi sites have a mirror plane, while Mn sites have inversion center as can be seen from Fig. 1. The valence atomic configuration of Bi is 6s 2 6p 3 and that of Mn is 3d 5 4s 2 .
Since strong correlations are suggested to play an important role in the properties of MnBi, we systematically study the effect of different exchange correlations (LDA and GGA) and U parameter values (0-4 eV) on the structure and magnetism. A rotationally invariant implementation of U is used in the calculations 16 and the results are summarized in Fig. 2.
We find that, in the U 5 0 eV limit, both LDA and GGA underestimate the volume with 16.6% and 4.4% errors respectively. The better agreement between GGA and experiments is because it gets the a lattice parameter more or less correct, while LDA underestimates both. Magnetic moments of Mn atoms are better with GGA and including a small U of 1-2 eV in GGA 1 U calculations improves c and m Mn , but agreement in a worsens. On the other hand, in the case of LDA, Hubbard U of more than 4 eV is necessary to reach good agreement with experiments, as can be seen from Fig. 2.
For comparison, we also calculated the optimized lattice parameters and magnetic moments using more accurate all-electron code WIEN2k within LDA, GGA-PBE and GGA-PBEsol approximations. The results are summarized in Table I. The LDA and GGA results are consistent with those form VASP calculations in Fig. 2. Interestingly, PBEsol does not substantially improve the results in this system. With spin polarization, the up and down d spin states are split by an exchange field of ,3.5 eV, as evident from the density of states in Fig. 4. The majority channel of Mn-d orbitals are fully occupied, however the minority channel is partially occupied by ,0.8 electrons     per Mn, leading to d 5.8 configurations at the Mn sites and a magnetic moment close to 4 m B . An opposite spin-polarization of about 20.15 m B is induced at the Bi sites as can be seen from the middle panel of Fig. 4. The calculated band structure and DOS compare reasonably well with previous calculations 10,14 . Also, experimental measurements using photoelectron spectroscopy 18 suggest that the occupied electron density peaks around 23 eV below the Fermi level, which is in agreement with Fig. 4. We find that, when spin-orbit coupling is included in the calculations, the energy of the system is lower by 1.2 meV/cell with the magnetic moments pointing along the a axis than when moments point along c axis. This leads to a MAE constant K u 5 21.97 MJ/m 3 , which is larger than the experimental value of 20.2 MJ/m 3 , but consistent with earlier theoretical calculations within GGA. Repeating the calculation for the moments pointing along [101] direction and using the expansion of anisotropy for hexagonal crystals, viz., K 5 K 1 sin 2 h 1 K 2 sin 4 h 1 …, we can estimate the leading terms. We get, K 1 5 21.2 meV and K 2 5 0.009 meV, which show that K 2 is too small to play any significant role in this system. Tight-binding model. A linear combination of atomic orbitals are used as the basis for constructing the tight-binding (TB) Hamiltonian and the periodic nature of the crystal is taken into account through the Bloch sum. From first-principles calculations discussed in the previous section, we know that the bands near Fermi level are made up of predominantly Bi-p and Mn-d states. Since there are two formula units in the primitive cell, the Hamiltonian spans 16 3 16 in orbital space and 2 3 2 in spin space. It has the following three component form, The first term, H K contains the electron kinetic energy and has same form for both spins. It can be written in terms of onsite energies e a and Slater-Koster hopping integrals E a,b (r ij ) within two-center approximation. Denoting the basis functions as jiasae, where i, a and s correspond to site, orbital and spin indices, we can write where, k is the wave vector and r ij is the vector connecting the neighbors i and j. For MnBi, we considered contributions from all neighbors and found that in addition to the nearest neighbor d 2 p coupling arising from Mn-Bi bonds, we also need to consider Mn-Mn interaction along c axis, nearest and next-nearest neighbor Bi-Bi interactions as shown in Fig. 5. Second neighbor coupling for Mn are weak and can be ignored to the lowest order. After performing the sum over the neighbors indicated in Fig. 5, the second part of Eq. 2 can be written as, The hopping matrices h dd , h dp and h pp are given in the Appendix. They depend on the wave vector k and the seven Slater-Koster parameters: t dds , t ddp , t ddd , t dps , t dpp , t pps , t ppp . The onsite energies e i are taken to be equal. Ignoring any spin-dependence of these parameters, we use ab-initio calculations with the paramagnetic structure to derive them. For hopping parameters, we assume a power-law distance dependence, t(r) 5 t(r 0 )(r 0 /r) n , where r 0 corresponds to bond lengths in the experimental structure. The exponent n can be estimated by fitting to DFT bands at different volumes. The values thus obtained are listed in Table II. The TB band structure calculated using the derived parameters are shown in Fig. 6 where the important features are found to agree well with DFT bands in Fig. 3. We can also see from Fig. 6(a), where the effect of Mn are ignored, that the large dispersion in Bi-p states arise from direct coupling between Bi atoms. The dispersion in d z 2 (flat line just above the Fermi level) is not well captured by the model, because we have ignored in-plane Mn coupling.
The magnetism is incorporated in the Hamiltonian via the second term in Eq. 1, which can be written in the Stoner form ignoring any k dependence, wherem is the direction of magnetization and t are the Pauli spin matrices. In this definition, the exchange splitting of each band is equal to the Stoner parameter I ia . In the present case, there are two parameters I d and I p , that correspond to splitting in d and p states respectively. They can be chosen to fit the experimental magnetic  With the parameters discussed above, the tight-binding model is able to reproduce the correct low temperature magnetic ground state. We get MAE of 22.7 meV, which is overestimated but has the correct sign. In solids, the strength of spin-orbit coupling is weakened by bonding 19 , which may be the source for this error. Thus, we show that a simple model involving only Bi-p and Mn-d orbitals can be used to study the anisotropy in these systems.
Lattice dynamics. Experimental measurements on MnBi have found that upon heating, the magnetic moments spontaneously change their direction from ab plane to c axis. Theoretical studies have suggested that this MAE anomaly can be explained by the lattice expansion upon heating 12,13 . The dynamics of the atomic lattice and the magnetic moments can also contribute towards the changes in MAE at high temperatures. To evaluate this, we investigate the lattice dynamics in MnBi in this section.
As shown in Fig. 7, the phonon band structure and DOS plots display a near complete separation of the Bi and Mn atom modes, with the Bi modes active between 0 and 2.5 THz and the much lighter Mn modes at higher frequencies between 4 and 6 THz. It is suggestive of the two systems being vibrationally decoupled for nearly all modes except the low frequency acoustic modes, which necessarily involve both atoms. It is of interest to examine this result in light of the mean squared displacements for the two atoms plotted in Figure 8. Despite the general separation of Mn and Bi modes, the mean squared displacements, for both atoms and both directions, are surprisingly consistent, varying by no more than 25 percent at the largest temperatures. This is in agreement with comparable force constants applicable to the Mn and Bi atoms, whose non-acoustic mode frequencies are then split due to the mass difference. Calculated mean square displacements U 11 and U 33 plotted in Fig. 8 show that the displacements of Bi atoms are similar in both a and c directions. However, the Mn atoms have higher displacements in the a direction. At room temperature, the mean amplitude of vibrations is around 0.1 Å .
In general, the high lying Mn optic modes have little dispersion, with the exception of one mode between C and A. The low energy modes are of more interest. Surprisingly, in the planar direction C to M, the two transverse modes are not nearly degenerate; in fact the    Continuing on the theme of elastic anisotropy, in Fig. 9 we present the calculated sound speeds for all three modes, as a function of propagation direction. These sound speeds were generated from first principles calculations of the elastic constants of MnBi using the all electron code WIEN2k, whose calculated elastic constants were found to be (in GPa) c 11 5 69.4, c 12 5 35.7, c 13 5 29.9, c 33 5 55.7, c 44 5 36.9. Several unusual features are found: one transverse mode, (a) has larger sound speeds in the axial direction than in plane, as is evident in Fig. 7. A second transverse mode (b) has a pointed structure along the axial direction, with the remainder of the surface essentially an ellipsoid. Finally, the longitudinal mode (c) appears to show its maximum sound speed in neither the planar nor c-axis direction.
Effect of lattice vibrations on MAE. At higher temperatures, all the phonon modes are excited and the lattice vibrations can be thought of as isotropic with mean displacements as shown in Fig. 8. Thus, it is possible to get a reasonable estimate of the magnetic anisotropy at different temperatures by averaging over configurations in which the atoms are displaced from their mean positions such that their mean square correspond to Fig. 8 at that temperature.
To make this problem computationally tractable, we consider only a single atomic displacement at a time and only along a and c directions and calculate the change in anisotropy energy K u as a function of deviations from the equilibrium positions, denoted as x. To minimize the effect of neighbors, we construct a larger 2 3 2 3 1 supercell. Calculated changes in total energies and MAE are shown in Fig. 10, where it can be seen that displacements upto 0.1 Å fall within thermal fluctuations, since kT , 25 meV. While all deviations increase the anisotropy, in-plane displacements have a stronger effect.
From Fig. 10, we see that DK u can be approximated as a quadratic function of x, so that K u (x) 5 K u (0) 1 kx 2 , where the coefficient k depends on the atom types and direction of x. Now, to obtain the average MAE, AEK u ae, at a particular temperature, we assume the atomic vibrations to be harmonic, i.e., x(t, T) 5 a(T) sin vt, where a(T) is the amplitude and v is the frequency of the vibrations. Then the average is defined as, where P(x) is the normalized probability function of a harmonic oscillator, Evaluating the integral yields, AEK u ae 5 kja(T)j 2 /2. Taking the ja(T)j 2 from Fig. 8 and summing over atoms and directions, we get AEK u ae as a function of temperature which is plotted in Fig. 11. Interestingly, it shows that atomic vibrations alone can result in spin-reorientation, albeit at a much higher temperature of 450 K. This simplified model shows that the change in MAE upon heating to 300 K, is ,0.8 meV, which is about 64% of the total change observed experimentally. Even though we have ignored all anharmonic effects and considered only isolated displacements, already these calculations suggest that lattice dynamics play an important role in the magnetic behavior of MnBi.

Discussion
First-principles electronic structure, as well as lattice dynamics calculations have been performed in an attempt to understand the magnetic properties of the ferromagnet MnBi. A tight binding model of this system has also been developed. Our calculations suggest unusual magneto elastic properties in this system. The calculated phonon spectra and sound speeds show anisotropic features, including a transverse mode with nearly double the planar sound speed of   the other transverse mode, as well as longitudinal sound speeds which are maximal at non-symmetry locations. All these observations together suggest that MnBi's magnetic interaction with the crystalline lattice is substantially more complex than previously believed. We also find that lattice vibrations contribute significantly to the temperature dependence of magnetic anisotropy and should be taken into account. This rich behavior including both magnetoelastic (bonding related) and magnetostrictive (anisotropy induced) couplings suggests a variety of interesting physical behavior such as temperature dependent magnon-phonon couplings, elastic and phonon anomalies under field and with temperature and potentially opportunities for tuning the magnetic behavior of MnBi. We believe experimental measurements of temperature dependent vibrational spectra are needed to gain further understanding of these issues. We note that, while the results imply that there could be interesting isotope effects on the magnetic behavior, Mn and Bi do not have convenient stable isotopes for such studies.

Methods
For first-principles calculations, we use DFT method implemented in the Vienna ab initio simulation package (VASP) 20,21 with a C-centered grid containing more than 11 k points in the Brillouin zone and 400 eV plane-wave energy cutoff. The gradient corrected GGA exchange-correlation functional is used for the calculations, unless otherwise noted 22 . Lattice parameters in Table I, elastic constants and sound velocities are calculated using the all-electron code WIEN2k 23 .
To obtain the phonon band structure shown in Fig. 7, we employed the PHONPY code 24 . It works by doing several total energy calculations on supercells incorporating ''frozen phonons,'' or atomic displacements dictated by the hexagonal crystal symmetry. Note that spin-orbit coupling is not incorporated in these calculations.
The sound velocities n in the solid are related to the elastic constants e and density r as, n~ffi ffiffiffiffiffiffi e=r p . The e is a 6 3 6 matrix with five independent parameters for hexagonal structures. The Christoffel equation is used to solve for the velocities.
In this paper, the magnetic anisotropy is defined as, K u~Eâ f g{Eĉ f g ð Þ where Eâ f g and Eĉ f g are the total energies of the system with magnetic moments pointing along the directions a and c respectively.