Magnetostructural coupling from competing magnetic and chemical bonding effects

An understanding of the origins of giant magnetostructural coupling is developed for the compound MnAs, a magnetic material that has served as a prototype for many promising technologies including caloric refrigeration, magnetic actuation, and spintronics. We demonstrate that strong coupling between magnetism and crystal structure arises from an orbital-speci(cid:12)c competition between exchange energies and kinetic (bonding) energies and that thermally activated spin-(cid:13)uctuations drive the unusual (cid:12)rst-order phase transition from high to low symmetry upon heating. The underlying mechanism raises the prospect of an exotic paramagnetic state featuring local (cid:13)uctuations in atomic positions and bonding on the time scale of the moment (cid:13)uctuations. The results should inform the design of new materials with enhanced magnetostructural coupling, found at the border between structural and magnetic stability.

Since the 1950s [1][2][3][4][5][6], MnAs has served as the prototypical material for strong coupling between magnetism and crystal structure. Upon heating, a dramatic firstorder magnetostructural phase transition is observed at T C = 318 K, whereupon the ferromagnetic (FM) hexagonal structure transforms to a lower-symmetry paramagnetic (PM) orthorhombic structure with a 2.5% smaller unit cell volume. As the temperature is further increased, the hexagonal structure is recovered via a continuous transition around T t = 398 K. The transition at T C is of great interest because it can be actuated by temperature, field, and mechanical stress. In particular, a moderate magnetic field can be used to drive giant entropy changes in MnAs, enabling magnetic refrigeration based on the magnetocaloric effect [7][8][9]. Other promising magnetocalorics, such as Gd 5 (Si,Ge) 4 , MnFe(P,X), and doped LaMnO 3 display similarly strong coupling [10][11][12][13][14][15][16]; however, the microscopic mechanism governing giant magnetostructural coupling in MnAs and related materials remains an open question.
In the ground state structure of MnAs, Mn is arranged on a hexagonal lattice [ Fig. 1(a)], with large spacing (3.72Å) within the basal plane and short spacing (2.85Å) perpendicular to the plane. Upon heating above 318 K or applying pressures above 2 kbar [17], first-order transitions to orthorhombic symmetry are observed [ Fig. 1(b)]. In these transitions, alternating rows of Mn move towards each other within the plane to create Mn-Mn zigzag chains with shorter contacts lengths (3.38Å at ambient pressure), while alternating rows of As atoms move a similar distance in the perpendicular direction. We refer to the magnitudes of these modes within the orthorhombic cell as δ Mn = ∆ Mn z /c and δ As = ∆As x /a, indicating deviations of the atoms' fractional coordinates from their ideal hexagonal positions. At 318 K and ambient pressure, δ Mn and δ As are 0.027 and 0.025, respectively [18].
There have been numerous attempts to explain the behavior of MnAs [4][5][6][19][20][21][22]. Goodenough and coworkers proposed that the transitions are Mn 3+ high-spin to low-spin transformations [6,17,23]. Conversely, Bean and Rodbell proposed a phenomenological model where a first-order transition can be realized in a sufficiently compressible material with a strongly volume-dependent magnetic exchange [4,5]. More recent computational efforts have concluded that the orthorhombic structure above T C is actually antiferromangetic (AFM), with the Mn moments alternating direction along the zigzag chains [C-type AFM, Fig. 1(c)] [19,[24][25][26]. However, neutron diffraction points to a PM state without long-range order [3], and paramagnetic scattering experiments [27,28] and chemical substitution studies [29] indicate ferromagnetic correlations and a local moment magnitude that remains approximately constant across T C . Furthermore, a variety of unusual physical properties reported in MnAs are largely unexplained, including the large amounts of local atomic disorder [30,31], anomalous elastic properties [32,33], and extremely low thermal conductivity [34]. Microscopic understanding consistent with the experimental situation has, so far, clearly been missing.
Here, we use electronic structure calculations on ordered and disordered magnetic states to establish physical mechanisms for both the pressure-driven and temperaturedriven phase transitions in MnAs. Both transitions are found to arise from orbital-specific competition between kinetic energy and different types of exchange energies: intra-atomic Hund's coupling for the pressure-driven transition, and interatomic magnetic exchange for the temperature-driven transition. Furthermore, our calculations suggest that this competition may become dynamic in the paramagnetic phase as magnetic fluctuations couple with large fluctuations in interatomic spacing. This is expected to result in an unusual disordered phase, which can explain the anomalous elastic and transport behavior.
We began our investigation by performing density functional theory-based calculations (DFT) [35][36][37] on MnAs using a generalized gradient approximation (PBEsol [39]). Starting with the experimental ambient pressure hexagonal [40] and orthorhombic [18] structures and ferromagnetic moments, we fully relaxed the lattice parameters and internal atomic coordinates. This resulted in cells with volumes respectively 13% and 19% smaller than experiment, in contrast with the typical slight overestimation of volumes often associated with the PBEsol approximation [41]. Additionally, the DFT energy of the orthorhombic cell was found to be lower than that of the hexagonal cell by 181 meV f.u. −1 , a result that is clearly unphysical given that the experimental ground state is ferromagnetic and hexagonal. Previous DFT studies have not addressed this inconsistency [25,42], presumably because a full relaxation into the global orthorhombic minimum was not investigated.
Employing the DFT+U electron correlation correction via Dudarev's U − J = U eff formalism [43] localizes the Mn d states that are systematically too diffuse in standard DFT [44] and brings both the lattice parameters and magnetic moments (3.4 µ B per Mn) of the hexagonal cell into good agreement with experiment for U eff = 1.2 eV with the PBEsol functional [ Fig. S1,S2]. Perhaps more importantly, this U eff proves to be the critical variable to correct the energetics, and the hexagonal cell is now 92 meV f.u. −1 lower in energy than the orthorhombic one. We also tried calculations with the standard PBE functional , and found very similar results: without U , the volumes are underestimated and the ground state is incorrect. The energetics and volume, however, can be corrected with a small U eff = 0.5 eV, which gives similar results to PBEsol with U eff = 1.2 eV [ Fig. S3].
Proceeding with PBEsol and U eff = 1.2 eV, we found that the metastable orthorhombic state has a smaller Mn moment than the hexagonal ground state (2.4 µ B vs. 3.4 µ B ). This moment and the DFT lattice parameters closely match (within 2 %) the experimental state obtained at low temperatures by applying pressure [45]. This indicates that the pressure-driven transition is consistent with a high-moment to low-moment transition. However, the volume of this orthorhombic cell is still 16% smaller than the ambient pressure orthorhombic cell observed at 318 K, and the δ Mn is much larger (0.061 vs. 0.027); therefore this low-moment cell can only explain the pressure-driven transition, not the thermally-driven transition.
We further investigated the energy surface of ferromagnetic MnAs during its transition from high-moment hexagonal to low-moment orthorhombic as a function of δ Mn and δ As . Figure 2(a) shows a 2D energy surface created from 900 calculations where Mn z and As x positions were held fixed, while all other coordinates and the lattice parameters were allowed to freely relax. We obtain a double-welled energy surface with local minima at the hexagonal structure (δ Mn = 0, δ As = 0) and in the low-moment orthorhombic structure (δ Mn = 0.061, δ As = 0.053), with the barrier between these wells indicated with a white dashed line running principally vertically on the phase diagram. This analysis shows that the transition is driven by the Mn sublattice: moving the Mn atoms without movement of the As atoms can traverse the phase transition line, but not vice-versa. We therefore proceed to plot the 1D energy surface and local Mn moment as a function of δ Mn only [ Fig. 2 where the As x position is now allowed to freely relax along with the other structural parameters. Here, a sharp drop in moment from 3.4 µ B to 2.4 µ B is observed as the hexagonal-to-orthorhombic phase boundary is crossed.
When we instead impose the C-type AFM structure, the energy surface looks very different. The potential is now single-welled, with an energy maximum at δ Mn = 0 and minimum at δ Mn = 0.05. Additionally, the local moment is virtually independent of δ Mn . The magnetic stabilization energy (FM minus AFM), which is proportional to Mn-Mn magnetic exchange energy along the zigzag chain, is at its most negative at δ Mn = 0, and changes sign to positive (favoring AFM) for δ Mn > 0.035, consistent with a previous extended Heisenberg model investigation [42]. Notably, the experimental ambient-pressure 318 K orthorhombic cell has δ Mn = 0.027, and δ Mn decreases as temperature rises. [46] Therefore, ambient-pressure MnAs always falls within the regime where ferromagnetic exchange dominates (δ Mn < 0.035), corroborating the experimental observation that the paramagnetic correlations are ferromagnetic, not antiferromagnetic. [27,28] At ambient pressure and 318 K, the observed δ Mn = 0.027 is close to the crossover in FM and AFM energies, consistent with the formation of a spin-fluctuationbased paramagnetic state. Since the dynamics of such a state are out of reach for DFT calculations, we simulate the paramagnetism using a special quasirandom structure (SQS) methodology [47], which has recently been applied successfully to paramagnets [48][49][50][51]. This approach involves constructing a supercell [we used the 96-atom supercell shown in Fig. 5(a)] and adjusting quasirandom up and down moments on the Mn atoms such that the local spin-spin correlations match the correlations of an infinite paramagnet (see Supplemental Material for details). When this disordered local moment cell is fully relaxed [gray star in Fig. 2 , it adopts an orthorhombic structure with the individual Mn and As atoms each relaxing differently in accordance with the broken local symmetries of the disordered cell. Overall, the average δ Mn is 0.035 ± 0.006. Even though this SQS calculation represents a static model for a dynamic phenomenon (paramagnetism), the agreement with experiment is remarkably good: the fully-relaxed SQS cell has an average local magnetic moment of 3.2 µ B per Mn, comparing well with 3.1 µ B from paramagnetic neutron scattering [27]. The relaxed lattice parameters are also all within 3% of their experimental 318 K values [18] (Table S2). Fig. 2(b-c) also shows an energy surface for the disordered state, calculated using a smaller 64-atom SQS by fixing the δ Mn values but allowing all other structural degrees of freedom to relax. This surface shows a very shallow minimum around δ Mn = 0.035 and a maximum at δ Mn = 0. We note that, by construction, the SQS calculations give higher energies than the ordered states, as the entropy which stabilizes the paramagnet is not considered in the calculated energy.
From these calculations, we can see that it is the high-moment ferromagnetism that establishes the highsymmetry lattice: any pressure or temperature sufficient to disrupt this magnetic configuration towards PM, AFM, or low-moment FM states will simultaneously drive a structural distortion to orthorhombic. To understand the origins of this behavior, we investigated the orbital bonding properties of MnAs. While Mn sits in a nearly perfect octahedron in the hexagonal structure, we found that the expected octahedral crystal field splitting is absent from the electronic structure [Supplemental Material Fig. S5]. This indicates that direct Mn-Mn interactions beyond the immediate Mn-As coordination shell play a significant role in the electronic structure, reminiscent of the FeAs family of superconductors which also exhibit strong and unusual magnetostructural coupling [52,53]. We therefore reoriented the d orbital basis according to the Mn-Mn contacts, as illustrated in Figs. 3(d) and Fig. 4(a). In this scheme, the d z 2 orbital lobes point out of the plane, along the nearest-neighbor Mn-Mn contact, while one lobe of the d x 2 −y 2 orbital points along the Mn-Mn zigzag contact that is activated by the dis- tortion. In the hexagonal structure, d x 2 −y 2 and d xy are degenerate, and d yz and d xz are degenerate. As shown in Fig. 3, all of the Mn d orbitals except d z 2 show partial DOS that are approximately half-filled with the majority spin states localized about 2.5 eV below the Fermi level with unoccupied minority spin states above. This electronic structure is stabilized by Hund's coupling (i.e. intra-atomic exchange), which is the driving force for the splitting of majority and minority spin states. Since only the majority spin states are occupied, this splitting lowers the energy. Conversely, the d z 2 orbital shows a much more dispersive DOS consistent with a band formed by Mn-Mn bonding along the short contacts Fig. 3(d). For this orbital, it is the kinetic energy benefit associated with chemical bonding that stabilizes the electronic structure.
In our projection scheme, the d x 2 −y 2 orbital shows large changes in projected moment when a Mn distortion is introduced (Fig. 4.), while the other d orbitals are comparatively unaffected (Fig. S6). In the ferromagnetic case, as the distortion forces Mn atoms together, overlap of the neighboring d x 2 −y 2 orbitals makes full spin-polarization untenable due to the Pauli exclusion principle, as electrons with the same spin cannot overlap spatially. This precipitates an electronic transition, seen in the reorganization of the DOS from a fully spin-polarized, localized structure in Fig. 4(b) to a much more dispersive structure in Fig. 4(d). The new electronic structure loses the Hund's coupling stabilization, but partially compensates by forming Mn-Mn bonds along the zigzag chain [ Fig. 4(e)]. The competition between Hund's coupling and bonding is clearly seen in Fig. 4(c), which shows the d x 2 −y 2 moment and the Mn-Mn zigzag contact bond strength (estimated using the integrated crystal orbital Hamilton population (-iCOHP) [54][55][56][57][58]) as a function of δ Mn . As δ Mn increases, the moment drops precipitously, while the -iCOHP rises from 0 eV to 0.6 eV, indicating the formation of a Mn-Mn bond. Furthermore, the orthorhombic distortion also increases the Mn-As and As-As -iCOHPs [ Fig. S7] as the asymmetric coordination optimizes bond distances. Altogether, these favorable bonding interactions explain the local stability of the low-moment orthorhombic state. Therefore, we can see that the high-moment hexagonal to low-moment orthorhombic transition experienced by MnAs under pressure can be understood as arising from a competition between Hund's coupling and chemical bonding. At ambient pressure, the Hund's coupling stabilization outweighs the kinetic energy consideration. However, the application of pressure forces a reorganization of the electrons to comply with the Pauli principle, destabilizing the high-volume high-moment state in favor of a lower-moment, more strongly bonded structure.
Even though the AFM state is never observed experimentally, it is informative to consider its behavior in order to understand how spin flips modify the system [Fig 4(f-i)]. In the AFM case, the Mn atoms that are brought together by the distortion have opposite spins and the Pauli exclusion principle does not come into play. As the d x 2 −y 2 orbitals begin to overlap, they are able to hybridize and form a metal-metal bond without disrupting the spin polarization. For this reason, the Mn local moment magnitudes remain approximately constant as a function of δ Mn , and bonding and moment can all be simultaneously optimized at δ Mn = 0.05. However, this AFM state is not the ground state due to non-optimized inter -atomic magnetic exchange [i.e. FM − AFM energy in Fig. 2(b)]. As discussed previously, this exchange is optimized at δ Mn = 0, stabilizing the FM hexagonal ground state.
Based on the analysis of the FM, AFM and SQS states, we are now equipped to understand why the ambientpressure transition from ferromagnetic to paramagnetic is coupled to a large structural distortion. The hexagonal ground state is not optimized from a bonding perspective, but maximizes both interatomic magnetic exchange energy and Hund's coupling and is therefore lower in energy than all competing states. As temperature rises and spin-flip fluctuations begin to set in, the magnetic exchange energy benefit is weakened, and, at 318 K, it becomes more favorable to structurally distort, forgoing the interatomic exchange energy in favor of the kinetic energy (bonding) stabilization of the orthorhombic state. The spin fluctuations aid this transition by introducing AFM pairs of Mn atoms along the zigzag chains, allowing d x 2 −y 2 orbitals to partially maintain their Hund's coupling stabilization when distorting. Interestingly, all the competition takes place in the d x 2 −y 2 orbital channel. The small Mn-Mn contact length along the d z 2 orbital causes it to remain in the bonded state at all temperatures, while long Mn-Mn contacts along the other orbitals cause them to stay in their spin-polarized Hund's coupled states at all temperatures. Therefore, the key to this transition is the competition between exchange and bonding energies in half-filled Mn-Mn orbitals that are just on the verge of possible overlap.
This interplay between magnetic exchange, magnetic moment, and structural distortion has interesting implications for the paramagnetic state. In the fully-relaxed SQS cell, we observed that the Mn and As atoms relaxed away from their symmetrical sites by up to 0.2Å, leading to very different local environments for different atoms. Such large local distortions have not been seen in other paramagnets simulated using the SQS method [48][49][50][51]. Figure 5(a) shows the relaxed SQS cell, and Fig. 5(b) shows all of the atoms from this supercell projected into a single primitive unit cell. Here, it can be seen that Mn and As atoms cluster into elongated clouds along the hexagonal-to-orthorhombic distortion modes; evidently, the PM state is most stable with a distribution of δ Mn and δ As values. In fact, this unusual finding can be connected to the mechanism of magnetostructural coupling described above.  than FM pairs, and bond strengths about twice as large. This is consistent with the idea that the Pauli exclusion principle disfavors d x 2 −y 2 orbital overlap in FM pairs of Mn atoms, but not in AFM pairs. Importantly, the other Mn-Mn contacts in the structure, along other d orbitals, show no dependence of the bond lengths or strengths on the spin correlation [ Fig. S8], highlighting once again that the competition-driven coupling occurs selectively in the Mn d x 2 −y 2 orbitals.
The dependence of Mn-Mn bonding on local spin configuration implies that in the paramagnetic state of MnAs, thermal spin fluctuations can be expected to couple to large atomic displacements accompanied by the dynamic formation and breaking of metal-metal bonds. We predict this exotic disordered state should be detectable by large and dynamic local-symmetry breaking [30,31], as well as strong magneto-vibrational inelastic neutron scattering, and anomalous thermal and spin transport. In particular, the local atomic disorder caused by this state may be invoked to explain the extremely low thermal conductivity [34] and the abnormal elastic properties [32,33]. Furthermore, the mechanisms for magnetostructural coupling driven by competition between exchange and kinetic energy scales described here are expected to be a general feature of other systems displaying strong coupling, and similar features of the disordered magnetic states may be expected.
This work was supported by the National Science Foundation through the DMR-1710638. Partial support to SDW through the NSF Materials Research Science and Electronic structure calculations were performed using density functional theory (DFT), as implemented in the Vienna Ab initio Simulation Package (VASP) [1] with projector augmented wave (PAW) pseudopotentials [2,3] within the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) revised for solids (PBEsol) [4,5]. A DFT+U approach was adopted, using the U − J = U eff formalism of Dudarev [6], with the U eff applied to the Mn d orbitals.
Once the appropriate U eff = 1.2 eV value was chosen (see below), calculations of the orthohexagonal unit cell [main text Fig. 1(b)] with varying amounts of Mn distortion and As distortion were performed. During the experimentally observed structural distortion from hexagonal to orthorhombic, the Mn atoms move principally along the c lattice direction, while the As atoms move along the a direction. Therefore, the parameters δ Mn = ∆Mn z /c and δ As = ∆As x /a were used to parameterize the distortion space. 900 calculations were performed on a 30 × 30 grid of δ Mn and δ As values ranging from -0.03 to 0.08. The Mn z parameter and As x parameter were kept fixed, while all of the other atomic coordinates and the unit cell were allowed to freely relax. Structural relaxations with an ionic convergence criterion of 1 × 10 −4 eV were performed three times iteratively to ensure convergence, and then static calculations with an electronic convergence criterion of 1 × 10 −6 eV were performed. Monkhorst-Pack k-point grids of 6 × 8 × 6 were used, and all calculations were initialized with ferromagnetic moments of 3 µ B on each Mn atom. Spin-orbit coupling was not included. In all calculations, the final configurations maintained P nma symmetry and ferromagnetic ordering.
Similarly, linear potential surfaces only as a function of δ Mn were calculated, with the As z coordinate now being allowed to relax freely. For these calculations, an increased k-point grid of 6 × 11 × 5 was used. The lattice parameters were initialized by linearly interpolating between the lattice parameters of the fully relaxed hexagonal structure (δ Mn = 0) and the local minimum energy orthorhombic structure, and then allowed to relax freely. These calculations were performed using both ferromagnetic and C-type antiferromagnetic moment initializations, and in all cases the magnetic ordering type stayed consistent with the initialization during relaxation. The results from each of these calculations were analyzed to obtain the crystal orbital Hamilton population (COHP) for nearby atoms and projected densities of states using the LOBSTER code [7][8][9][10][11]. The projected densities of states were calculated using a non-standard orientation of the d orbital basis set, as described in the main text and below. In order to do this, the static VASP calculations that fed into the LOBSTER calculations were performed with the MnAs structure rotated in 3d space such that correct Mn-Mn contacts were oriented along the cartesian z and x directions. pymatgen [12] was used to perform some of the cell manipulations and analysis, and VESTA [13] was used to visualize the structures and charge densities.

B. Simulation of paramagnetic state using special quasirandom structures
The paramagnetic state of MnAs was simulated using special quasi random structures (SQS) [14] obtained using mcsqs program [15] included within the ATAT package. A paramagnetic SQS is a supercell which is decorated with up and down moments on the Mn S2 atoms in such a way as to accurately reproduce the local spin-spin correlations expected in a random, dynamic paramagnet using just a single cell. A 96-atom SQS with supercell size 2 × 3 × 2 of the the primitive orthorhombic cell was generated such that the 12 smallest pair correlations (radius up to 6Å) and its first triplet correlation (radius 3.4Å) match those of the random state. A full structural optimization was performed on this cell, starting with lattice parameters and atom positions from the lowest-energy AFM cell. The Mn atoms were each initialized with magnetic moments of positive or negative 3 µ B based on the constructed SQS, and symmetry was switched off such that the individual Mn and As atoms were each allowed to move in any direction and locally break the P nma structural symmetry, and the cell was allowed to change shape arbitrarily. A Γ-centered k-point grid of 3 × 3 × 3 was used. This structural optimization was performed several times iteratively, always resetting the initialized magnetic moments, until a final force convergence of -0.005 meVÅ −1 for all atoms was reached. The atoms all moved substantially away from their P nma positions, and the quasirandom magnetic moments on each Mn atom remained stable with an overall cell moment of nearly zero. A static calculation with energy convergence of 1 × 10 −6 eV was then performed, and the COHP was calculated. The resulting cell is shown in table S1.
Additionally, a 64-atom 2 × 2 × 2 SQS was generated such that the 17 smallest pair correlations (radius up to 6.8Å) match the random state. This cell was used to approximately obtain an energy surface of the random state as a function of δ Mn by performing selective dynamics structural relaxations on the SQS with fixed values of the Mn z parameters while all other atomic coordinates were allowed to fully relax. Structural optimizations were performed using a Γ-centered k-point grid of size 2 × 4 × 2 and an energy convergence for the ionic loop of 0.001 eV. For each calculation in the energy surface, these optimizations were performed three times iteratively to ensure convergence.  S1. Dependence of the MnAs distortion energy landscape on the U eff = U − J parameter added to the Mn d-orbitals. With U eff = 0, the distorted orthorhombic low-moment structure is erroneously lower in energy than the true (experimentally observed) hexagonal ground state (δ Mn = 0). Upon adding the U correction, the hexagonal state is stabilized relative to the orthorhombic state and becomes lower in energy for U eff = 1.0 eV and above. The curve for the U value selected for use in our study (U eff = 1.2 eV) is shown in bold. It should be noted that the chosen value of U eff is dependent on the details of DFT calculations. For example, when we attempted to use the standard PBE GGA instead of PBEsol, the qualitative behavior of the energy curves remained the same but the optimal U eff was found to be 0.5 eV.    Fig. 4. However, since there are many more Mn-As and As-As contacts within the unit cell, the total changes in Mn-As and As-As -iCOHPs per cell are generally larger than the Mn-Mn contribution. From these calculations, it can be see that the overall bonding is optimized around δ = 0.06 for the FM case, and δ = 0.05 for the AFM case.

S11
VII. Additional bond-lengths and -iCOHPs in the relaxed SQS cell Unlike the zigzag contact, the two contacts shown here do not exhibit a dependence of bond length or -iCOHP on the spin correlation. This observation highlights that the competition-driven magnetostructural coupling is only active in a single Mn d orbital, the d x 2 −y 2 orbital discussed in the main text.