Magnon-polarons in cubic collinear Antiferromagnets

TL;DR

This paper provides a theoretical framework for understanding magnon-polaron excitations in cubic antiferromagnets, demonstrating how hybridization occurs and can be controlled, with a detailed case study on Nickel Oxide.

Contribution

It introduces a general approach to identify hybridizing magnon and phonon modes and applies it to Nickel Oxide, revealing unique features of antiferromagnetic magnon-polarons.

Findings

Hybridization between magnons and longitudinal phonons in antiferromagnets.

External magnetic fields can control magnon-polaron hybridization.

Qualitative analysis matches quantitative results for Nickel Oxide.

Abstract

We present a theoretical study of excitations formed by hybridization between magnons and phonons - magnon-polarons - in antiferromagnets. We first outline a general approach to determining which magnon and phonon modes can and cannot hybridize in a system thereby addressing the qualitative questions concerning magnon-polaron formation. As a specific and experimentally relevant case, we study Nickel Oxide quantitatively and find perfect agreement with the qualitative analysis, thereby highlighting the strength of the former. We find that there are two distinct features of antiferromagnetic magnon-polarons which differ from the ferromagnetic ones. First, hybridization between magnons and the longitudinal phonon modes is expected in many cubic antiferromagnetic structures. Second, we find that the very existence of certain hybridizations can be controlled via an external magnetic field, an effect which comes in addition to the ability to move the magnon modes relative to the phonons modes.

Figures

Equations

$$\begin{array}{rl}{\mathcal{H}}_{\mathrm{ME}}^{\mathrm{AFM}}& =\sum _{\alpha \beta }\int {\mathrm{d}}^{3}r{B}^{\alpha \beta }{n}^{\alpha }(\mathbf{r})n^{\beta }(\mathbf{r}){ϵ}^{\alpha \beta }(\mathbf{r})\\ & +\sum _{\alpha \beta }\int {\mathrm{d}}^{3}r{\tilde{B}}^{\alpha \beta }{n}^{\alpha }(\mathbf{r})n^{\beta }(\mathbf{r}){\widetilde{ϵ}}^{\alpha \beta }(\mathbf{r}),\end{array}$$

$${\mathcal{H}}_{\mathrm{ME}}^{\mathrm{circ}}=\frac{B^{\perp }{n}^z}{4}\left[\alpha \left(\frac{\partial u^x}{\partial r^z}-i\frac{\partial u^y}{\partial r^z}\right)+\beta \left(\frac{\partial u^x}{\partial r^z}+i\frac{\partial u^y}{\partial r^z}\right)\right].$$

$$\begin{array}{rl}{\mathcal{H}}_{\mathrm{ME}}^{\mathrm{NiO}}& =\frac{n^{\prime x}{n}^{\prime z}}{3\sqrt{2}}\left[2B^{\parallel }\frac{\partial u^z}{\partial r^z}+\frac{B^{\perp }}{2}\left(\frac{\partial u^x}{\partial r^z}+\frac{\partial u^y}{\partial r^z}\right)\right]\\ & +\frac{n^{\prime y}{n}^{\prime z}}{2\sqrt{3}}{B}^{\perp }\left(\frac{\partial u^y}{\partial r^z}-\frac{\partial u^x}{\partial r^z}\right),\end{array}$$

$$\begin{array}{rl}{\mathcal{H}}_{\mathrm{ME}}^{\mathrm{NiO}}& =\{\left(A\alpha +B\beta \right)\left(\frac{1}{3\sqrt{2}}\right)\left[2B^{\parallel }\frac{\partial u^z}{\partial r^z}+\frac{B^{\perp }}{2}\left(\frac{\partial u^x}{\partial r^z}+\frac{\partial u^y}{\partial r^z}\right)\right]\\ & +\left(A\beta -B\alpha \right)\left[\frac{i}{2\sqrt{3}}{B}^{\perp }\left(\frac{\partial u^y}{\partial r^z}-\frac{\partial u^x}{\partial r^z}\right)]\right\}\times \frac{n^{\prime z}}{(A^2+B^2)}.\end{array}$$

$$\begin{array}{r}\widetilde{\mathbf{ϵ}}=\left(\begin{array}{ccc}{ϵ}^{xz}& {ϵ}^{xz}+{ϵ}^{yz}& {ϵ}^{zz}\\ {ϵ}^{xz}+{ϵ}^{yz}& {ϵ}^{yz}& {ϵ}^{zz}\\ {ϵ}^{zz}& {ϵ}^{zz}& 0\end{array}\right).\end{array}$$

$$\begin{array}{rl}{\mathcal{H}}_{\mathrm{ME}}^{\mathrm{NiO}}& =\frac{\sqrt{2}{n}^{\prime x}{n}^{\prime z}}{3}\left[2{\tilde{B}}^{\perp }\frac{\partial u^z}{\partial r^z}-\frac{({\tilde{B}}^{\parallel }+{\tilde{B}}^{\perp })}{2}\left(\frac{\partial u^x}{\partial r^z}+\frac{\partial u^y}{\partial r^z}\right)\right]\\ & +\frac{n^{\prime y}{n}^{\prime z}}{\sqrt{3}}{\tilde{B}}^{\parallel }\left(\frac{\partial u^x}{\partial r^z}-\frac{\partial u^y}{\partial r^z}\right).\end{array}$$

$$\begin{array}{rl}{\mathcal{H}}_{\mathrm{m}}& =\frac{1}{2}\sum _{i,{\mathbf{\delta }}_1}{J}_1{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+{\mathbf{\delta }}_1}+\frac{1}{2}\sum _{i,{\mathbf{\delta }}_2}{J}_2{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+{\mathbf{\delta }}_2}\\ & +\sum _i{D}_{x^{\prime }}{\left(S_i^{x^{\prime }}\right)}^2+\sum _i{D}_{y^{\prime }}{\left(S_i^{y^{\prime }}\right)}^2,\end{array}$$

$$\begin{array}{rl}{S}_i^{\mathrm{A}{z}^{\prime }}& =S-a_i^{†}{a}_i^{},\\ S_i^{\mathrm{A}+}& =\sqrt{2S}{a}_i^{},\\ S_i^{\mathrm{A}-}& =\sqrt{2S}{a}_i^{†},\end{array}\begin{array}{rl}{S}_j^{\mathrm{B}{z}^{\prime }}& =-S+b_j^{†}{b}_j^{},\\ S_j^{\mathrm{B}+}& =\sqrt{2S}{b}_j^{†},\\ S_j^{\mathrm{B}-}& =\sqrt{2S}{b}_j^{},\end{array}$$

$$\begin{array}{rl}{a}_i^{}& =\frac{1}{\sqrt{N_A}}\sum _{\mathbf{k}}{e}^{-i\mathbf{k}\cdot {\mathbf{x}}_{\mathbf{i}}}{a}_{\mathbf{k}}^{},\\ a_i^{†}& =\frac{1}{\sqrt{N_A}}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot {\mathbf{x}}_{\mathbf{i}}}{a}_{\mathbf{k}}^{†},\end{array}\begin{array}{rl}{b}_j^{}& =\frac{1}{\sqrt{N_B}}\sum _{\mathbf{k}}{e}^{-i\mathbf{k}\cdot {\mathbf{x}}_{\mathbf{j}}}{b}_{\mathbf{k}}^{},\\ b_j^{†}& =\frac{1}{\sqrt{N_B}}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot {\mathbf{x}}_{\mathbf{j}}}{b}_{\mathbf{k}}^{†},\end{array}$$

$${\gamma }_{n\mathbf{k}}^{ab}=\sum _{{\mathbf{\delta }}_n\in ab}{e}^{i\mathbf{k}\cdot {\mathbf{\delta }}_n}.$$

$$\displaystyle\begin{split}\mathcal{H}_{\rm m}=&\sum_{\bm{k}}\Bigg{[}A_{\bm{k}}\left(a_{\bm{k}}^{\dagger}a_{\bm{k}}^{{\mathstrut}}+b_{\bm{k}}^{\dagger}b_{\bm{k}}^{{\mathstrut}}\right)+B_{\bm{k}}\left(a_{\bm{k}}^{{\mathstrut}}b_{-\bm{k}}^{{\mathstrut}}+a_{\bm{k}}^{\dagger}b_{-\bm{k}}^{\dagger}\right)\\ &\quad+C\left(a_{\bm{k}}^{{\mathstrut}}a_{-\bm{k}}^{{\mathstrut}}+b_{\bm{k}}^{{\mathstrut}}b_{-\bm{k}}^{{\mathstrut}}\right)+C\left(a_{\bm{k}}^{\dagger}a_{-\bm{k}}^{\dagger}+b_{\bm{k}}^{\dagger}b_{-\bm{k}}^{\dagger}\right)\Bigg{]},\end{split}$$

$$A_{\mathbf{k}}$$

$$B_{\mathbf{k}}$$

$$C$$

$$\begin{array}{r}{\mathbf{\xi }}_{\mathbf{k}}^{}=\left(\begin{array}{c}{a}_{\mathbf{k}}^{}\\ b_{-\mathbf{k}}^{†}\\ a_{-\mathbf{k}}^{†}\\ b_{\mathbf{k}}^{}\end{array}\right),\end{array}\begin{array}{r}{\mathbf{H}}_{\mathbf{k}}^{\mathrm{m}}=\left(\begin{array}{cccc}{A}_{\mathbf{k}}& B_{\mathbf{k}}& 2C& 0\\ B_{\mathbf{k}}& A_{\mathbf{k}}& 0& 2C\\ 2C& 0& A_{\mathbf{k}}& B_{\mathbf{k}}\\ 0& 2C& B_{\mathbf{k}}& A_{\mathbf{k}}\end{array}\right),\end{array}$$

$${\mathcal{H}}_{\mathrm{m}}=\frac{1}{2}\sum _{\mathbf{k}}{\xi }_{\mathbf{k}}^{†}{\mathbf{H}}_{\mathbf{k}}^{\mathrm{m}}{\xi }_{\mathbf{k}}^{}.$$

$$\mathbf{g}{\mathbf{H}}_{\mathbf{k}}^{\mathrm{m}}{\mathbf{T}}_{\mathbf{k}i}=\hslash {\omega }_{\mathbf{k}}{g}_{ii}{\mathbf{T}}_{\mathbf{k}i},$$

$$\hslash {\omega }_{\mathbf{k}}^{\pm }=\frac{1}{2}\sqrt{A_{\mathbf{k}}^2-{\left(B_{\mathbf{k}}^{}\pm 2C\right)}^2},$$

$$-{\omega }^2{m}_s{u}_{is}^{\alpha }=\sum _j\sum _t\sum _{\beta }{K}_{is,jt}^{\alpha \beta }{u}_{jt}^{\beta }.$$

$${\mathcal{H}}_{\mathrm{p}}=\sum _{\mathbf{q}\lambda }\hslash {\omega }_{\mathbf{q}\lambda }{c}_{\mathbf{q}\lambda }^{†}{c}_{\mathbf{q}\lambda }^{},$$

$$\begin{array}{rl}{\mathbf{ϵ}}_{\mathbf{q}1}& =[\cos {\theta }_{\mathbf{q}}\cos {\varphi }_{\mathbf{q}},\cos {\theta }_{\mathbf{q}}\sin {\varphi }_{\mathbf{q}},-\sin {\theta }_{\mathbf{q}}],\\ {\mathbf{ϵ}}_{\mathbf{q}2}& =i[-\sin {\varphi }_{\mathbf{q}},\cos {\varphi }_{\mathbf{q}},0],\\ {\mathbf{ϵ}}_{\mathbf{q}3}& =i[\sin {\theta }_{\mathbf{q}}\cos {\varphi }_{\mathbf{q}},\sin {\theta }_{\mathbf{q}}\sin {\varphi }_{\mathbf{q}},\cos {\theta }_{\mathbf{q}}],\end{array}$$

$${\mathcal{H}}_{\mathrm{ME}}=\sum _{\alpha \beta \gamma \lambda }\sum _{i,\mathbf{\delta }}{B}_{i,i+\mathbf{\delta }}^{\alpha \beta \gamma \lambda }{S}_i^{\alpha }{S}_{i+\mathbf{\delta }}^{\beta }{R}_{i,i+\mathbf{\delta }}^{\gamma \delta },$$

$$R^{\gamma\lambda}_{ij}=\frac{1}{2}\frac{1}{|\bm{r}_{i}-\bm{r}_{j}|^{2}}\Bigg{[}\left(r_{i}^{\gamma}-r_{j}^{\gamma}\right)\left(u_{i}^{\lambda}-u_{j}^{\lambda}\right)+\left(r_{i}^{\lambda}-r_{j}^{\lambda}\right)\left(u_{i}^{\gamma}-u_{j}^{\gamma}\right)\Bigg{]}.$$

$${\mathcal{H}}_{\mathrm{ME}}=\sum _{\alpha \beta }\sum _{i,\mathbf{\delta }}{B}_{|\mathbf{\delta }|}^{\alpha \beta }{S}_i^{\alpha }{S}_{i+\mathbf{\delta }}^{\beta }{R}_{i,i+\mathbf{\delta }}^{\alpha \beta },$$

$${\mathbf{u}}_i=\sum _{\mathbf{q},\lambda }{\mathbf{ϵ}}_{\mathbf{q}\lambda }\sqrt{\frac{\hslash }{2m{\omega }_{\mathbf{q}\lambda }N}}\left(c_{\mathbf{q},\lambda }^{†}+c_{-\mathbf{q},\lambda }\right)e^{i\mathbf{q}\cdot {\mathbf{r}}_i},$$

$$\begin{array}{rl}{R}_{i,i+\mathbf{\delta }}^{\alpha \beta }& =\sum _{\mathbf{q},\lambda }\sqrt{\frac{\hslash }{2m{\omega }_{\mathbf{q}\lambda }N}}\left(c_{\mathbf{q},\lambda }^{†}+c_{-\mathbf{q},\lambda }\right)e^{i\mathbf{q}\cdot {\mathbf{r}}_i}\\ & \times \left({\delta }_{}^{\alpha }{\mathbf{ϵ}}_{\mathbf{q}\lambda }^{\beta }+{\delta }_{}^{\beta }{\mathbf{ϵ}}_{\mathbf{q}\lambda }^{\alpha }\right)\left(1-e^{i\mathbf{q}\cdot \mathbf{\delta }}\right).\end{array}$$

$$G_{\mathbf{q}\lambda \mathbf{\delta }}^{\alpha \beta }=\frac{B_{\text{absolutevalue}\mathbf{\delta }}^{\alpha \beta }}{2{\mathbf{\delta }}^2}\sqrt{\frac{\hslash }{2m{\omega }_{\mathbf{q}\lambda }N}}\left({\delta }^{\alpha }{\hat{e}}_{\mathbf{q}\lambda }^{\beta }+{\delta }^{\beta }{\hat{e}}_{\mathbf{q}\lambda }^{\alpha }\right)\left(1-e^{i\mathbf{q}\cdot \mathbf{\delta }}\right),$$

$${\mathcal{H}}_{\mathrm{ME}}=\sum _{i,\mathbf{\delta }}\sum _{\alpha \beta }\sum _{\mathbf{q}\lambda }{G}_{\mathbf{q}\lambda \mathbf{\delta }}^{\alpha \beta }{S}_i^{\alpha }{S}_{i+\mathbf{\delta }}^{\beta }\left(c_{\mathbf{q}\lambda }^{†}+c_{-\mathbf{q}\lambda }^{}\right)e^{i\mathbf{q}\cdot {\mathbf{r}}_i}.$$

$${\mathcal{H}}_{\mathrm{ME}}=\sum _{i,\mathbf{\delta }}\sum _{\alpha \beta }\sum _{\mathbf{q}\lambda }{\tilde{G}}_{\mathbf{q}\lambda \mathbf{\delta }}^{\alpha \beta }{S}_i^{\prime \alpha }{S}_{i+\mathbf{\delta }}^{\prime \beta }\left(c_{\mathbf{q}\lambda }^{†}+c_{-\mathbf{q}\lambda }^{}\right)e^{i\mathbf{q}\cdot {\mathbf{r}}_i}.$$

$$\displaystyle\begin{split}\mathcal{H}_{\rm ME}&=\sqrt{\frac{N_{\rm A}S^{3}}{2}}\sum_{\bm{q}\lambda}\left(c^{\dagger}_{\bm{q}\lambda}+c^{\mathstrut}_{-\bm{q}\lambda}\right)\times\\ &\quad\Bigg{\{}\left(\sum_{\bm{\delta}\in\mathrm{AA}}\tilde{G}_{{\bm{q}\lambda}\bm{\delta}}^{xz}-\sum_{\bm{\delta}\in\mathrm{AB}}\tilde{G}_{{\bm{q}\lambda}\bm{\delta}}^{xz}\right)\left(a_{\bm{q}}^{\mathstrut}+a_{-\bm{q}}^{\dagger}-b_{\bm{q}}^{\mathstrut}-b_{-\bm{q}}^{\dagger}\right)\\ &\>-i\left(\sum_{\bm{\delta}\in\mathrm{AA}}\tilde{G}_{{\bm{q}\lambda}\bm{\delta}}^{yz}-\sum_{\bm{\delta}\in\mathrm{AB}}\tilde{G}_{{\bm{q}\lambda}\bm{\delta}}^{yz}\right)\left(a_{\bm{q}}^{\mathstrut}-a_{-\bm{q}}^{\dagger}+b_{\bm{q}}^{\mathstrut}-b_{-\bm{q}}^{\dagger}\right)\Bigg{\}}.\end{split}$$

$${\tilde{M}}_{\mathbf{q}\lambda }=\sqrt{\frac{N_{\mathrm{A}}{S}^3}{2}}\left[\sum _{\mathbf{\delta }\in \mathrm{AA}}\left({\tilde{G}}_{\mathbf{q}\lambda \mathbf{\delta }}^{xz}-i{\tilde{G}}_{\mathbf{q}\lambda \mathbf{\delta }}^{yz}\right)-\sum _{\mathbf{\delta }\in \mathrm{AB}}\left({\tilde{G}}_{\mathbf{q}\lambda \mathbf{\delta }}^{xz}-i{\tilde{G}}_{\mathbf{q}\lambda \mathbf{\delta }}^{yz}\right)\right],$$

Full text

I Introduction

Ultrafast magnetization dynamics is expected to play a pivotal role in the development of Terahertz (THz) technologies Walowski and Münzenberg (2016); Mittleman (2017). These processes cover the elusive THz frequency gap that roughly spans from 100 GHz to 30 THz Dhillon et al. (2017); Tonouchi (2007), the upper and lower limits of microwave and optical techniques, respectively. An early attempt on ultrafast magnetization dynamics was made in ferromagnetic Nickel almost 20 years ago using THz lasers Beaurepaire et al. (1996). Several new findings have propelled the field Kirilyuk et al. (2010); Koopmans et al. (2005), such as ultrafast coherent control of spin waves Kampfrath et al. (2011), ultrafast optically induced magnetization Ju et al. (2004); Battiato et al. (2010); Berritta et al. (2016), magnetic switching Stanciu et al. (2007); Radu et al. (2011); Stupakiewicz et al. (2017) and heat-assisted magnetic recording Sander et al. (2017). The complex interplay among charge, spin, elastic and optical degrees of freedom underlies the rich physics governing the ultrafast magnetization dynamics.

Antiferromagnetic materials provide a natural niche in this field Ivanov (2014); Jungwirth et al. (2016); Baltz et al. (2018); Jungwirth et al. (2018); Gomonay et al. (2018). Their fast magnetization dynamics, with the potential to cover the Terahertz range Baltz et al. (2018); Dhillon et al. (2017), and the lack of net magnetic moment Baltz et al. (2018), have instigated a growing interest in antiferromagnets (AFMs). AFM insulators are particularly interesting due to the absence of Joule heating caused by the scattering of charge currents Lin et al. (2016); Wu et al. (2016); Bender et al. (2017); Lebrun et al. (2018). The focus is rather on the spin currents carried by magnons, the quantized excitations of the magnetization dynamics. The ultrafast magnetization dynamics in AFM insulators can thus be understood in terms of these magnons and their interaction with the phonons - the quantized excitations of the lattice vibrations.

Although the exchange interaction underlies magnetism, the magnon-phonon interaction is crucial for the dynamics and equilibration of the spin system Akhiezer et al. (1968). The latter is brought on by the non-linear processes in which magnons scatter while absorbing or emitting phonons Akhiezer et al. (1968). In contrast, the linear magnetoelastic coupling Kittel (1949); Akhiezer et al. (1968) results in magnons and phonons combining to form hybrid excitations - magnon-polarons Kamra and Bauer (2014); Kamra et al. (2015) - when their coupled dispersions anticross. Despite the wavevector range corresponding to significant hybridization being small, it has been found to act as an energy short-circuit between the magnon and phonon subsystems Rückriegel et al. (2014); Bozhko et al. (2017). While the magnon-polarons have been studied in great detail Kittel (1949, 1958); Akhiezer et al. (1968), there has been a rekindling of interest in the phenomenon due to recent advances in fabrication and measurement techniques as well as fresh breakthroughs in the field of spintronics Chumak et al. (2015); Uchida et al. (2011); Saitoh et al. (2006). Despite decades of study, key questions like which phonon mode should or should not hybridize with magnons remain insufficiently understood. A related question is the issue of spin conservation in magnets Holanda et al. (2018), which is often invoked to understand several phenomena including the formation of magnon-polarons. In particular, it has been shown that spin conservation may not be invoked in a simple manner while addressing the effects of magnetoelastic interaction Kamra et al. (2015). This is because the latter is primarily rooted in spin-nonconserving interactions such as spin-orbit and dipolar contributions Kittel (1949); Akhiezer et al. (1968).

Magnon-polarons formed in the ferrimagnetic insulator yttrium iron garnet (YIG) have revealed their direct signatures in several recent experiments Flebus et al. (2017); Rückriegel et al. (2014); Man et al. (2017). Pronounced features, attributed to magnon-polarons, in the magnetic field dependence of the spin Seebeck effect (SSE) Kikkawa et al. (2016); Flebus et al. (2017) as well as nonlocal spin transport Cornelissen et al. (2017) have been observed in YIG. The groundwork for these observations was laid out in previous works. For instance, coherent elastic waves were used in spin pumping Weiler et al. (2012); Uchida et al. (2011) experiments, as well as spin wave excitations Kamra et al. (2015); Li et al. (2017). AFMs represent a step forward in this field. Although the ultrafast-response of AFM insulators is an exciting property, its control with magnetic fields is challenging. Thus its manipulation via magnetoelastic effects presents a useful alternative.

Here, we present a theoretical study of coupled magnetoelastic modes in AFMs. As compared to their ferromagnetic counterparts, the multisublattice nature of AFMs hosts qualitatively new, subtle, and rich magnetoelastic phenomena. Focusing on magnon-polaron formation, we outline a general method for gaining qualitative, physical insight into which magnon and phonon modes hybridize, given a crystal symmetry and ground state. As a special case, we study the widely used Nickel Oxide (NiO) in its collinear ground state, in which its spins are oriented at an angle with the crystal axes. This misalignment is found to permit novel effects such as hybridization with longitudinal phonons propagating along a crystal axis. Furthermore, NiO is found to host linearly polarized magnons which attain an increasingly elliptical polarization on application of an external magnetic field. This control of the qualitative nature of the magnons permits a magnetic field control of whether or not the magnons participate in forming magnon-polarons. This tunability goes well beyond the Zeeman energy shift afforded by the ferromagnetic magnons Kikkawa et al. (2016); Flebus et al. (2017), and opens prospects for novel functionalities.

This paper is structured as follows. In Section II, we present a general analysis of magnon-polarons in AFMs determining which phonons do or do not hybridize with which magnons. Considering NiO as an apt example, we examine its eigenmodes quantitatively finding them to be consistent with the general analysis. In Sec. III, we derive the quantum theory describing magnons, phonons and the magnon-phonon interaction. The quantitative results of the magnon-polarons states in NiO follows in Sec. IV, which in turn are compared with qualitative results predicted in Sec. II. Finally, we end with discussions and conclusions in Sec. V and VI, respectively.

II Semi-classical qualitative analysis

In the continuum (long-wavelength) limit, the magnetoelastic Hamiltonian in cubic AFMs is given as

[TABLE]

where $\epsilon^{\alpha\beta}=\frac{1}{2}\left(\frac{\partial u^{\alpha}}{\partial r^{\beta}}+\frac{\partial u^{\beta}}{\partial r^{\alpha}}\right)$ is the strain tensor, where $u^{\alpha}$ is the atom displacement field in the $\alpha$ direction. ${B^{\alpha\beta}=B^{\parallel}\delta^{\alpha\beta}+B^{\perp}\left(1-\delta^{\alpha\beta}\right)}$ and ${\tilde{B}^{\alpha\beta}=\tilde{B}^{\parallel}\delta^{\alpha\beta}+\tilde{B}^{\perp}\left(1-\delta^{\alpha\beta}\right)}$ are 4 magnetoelastic coefficients, and $\bm{n}$ is the Néel field. $\tilde{\epsilon}$ is an elastic tensor with elements being linear combinations of elements of the strain tensor. The first term is the antiferromagnetic analogy of the conventional ferromagnetic magnetoelastic Hamiltonian Kittel (1949). The second term, where $\tilde{\epsilon}$ appears, derives from the internal spin structure. This term is unique to materials with at least two sublattices, and the exact form of $\tilde{\epsilon}$ depends upon the spin structure of the AFM. If the spin structure is trivial, meaning that all $n$’th nearest neighbors of a lattice site belong to a single sublattice for any $n$, then this term disappears. We will refer to spin structures which do not fulfill this requirement as being non-trivial. A full derivation of this Hamiltonian is given in Appendix C.

To start with we will consider a simple AFM with a trivial spin structure, where only the first term in the antiferromagnetic magnetoelastic Hamiltonian (1) appears. We will only consider magnons and phonons propagating along one of the crystal axes, which we define to be the $\bm{\hat{z}}$ direction for concreteness. The three independent phonon modes are then proportional to $\partial u^{\gamma}/\partial r^{z}$, where $\gamma\in\{x,y\}$ describe transverse phonons, whereas $\gamma=z$ describes a longitudinal phonon.

In order to give a physical interpretation of the hybridization, we need to express the Hamiltonian in terms of both the free magnon and phonon eigenmodes. The strain tensor components $\epsilon^{\alpha\beta}$ are superpositions of phonon eigenmodes, and the Hamiltonian is thus already given in terms of phonon eigenmodes. The remaining task is therefore to find the magnon eigenmodes expressed as function of the Néel field.

II.0.1 Circularly polarized magnons in easy axis antiferromagnets

Let us now consider an easy axis AFM where the spins align along the $\bm{\hat{z}}$ axis in the (classical) ground state. In easy-axis AFMs the Hamiltonian is invariant to a global spin rotation about the $\bm{\hat{z}}$ axis, and we therefore expect the magnons modes to be circularly polarized. Moreover, as magnons corresponds to small deviations from the spin ground state, magnons leave $n^{z}$ approximately constant, whereas $n^{x}$ and $n^{y}$ are expected to oscillate harmonically. These considerations combined imply that we may express the magnon modes as $\alpha\equiv n_{x}+in_{y}$ and $\beta\equiv n_{x}-in_{y}$. The magnetoelastic Hamiltonian can then be rewritten as

[TABLE]

Note that in this situation, angular momentum in the $\bm{\hat{z}}$ direction is conserved. If we create a magnon, this may hybridize and produce a circularly polarized phonon, which has got angular momentum along $\bm{\hat{z}}$. If we create a linearly polarized phonon (with zero angular momentum along $\bm{\hat{z}}$), it may hybridize and produce spinless combination of magnons $\propto\alpha\pm\beta$. Angular momentum is hence conserved. This follows directly from the rotational symmetry about the $\bm{\hat{z}}$ axis 111Angular momentum conservation requires a continuous rotational symmetry about the $\bm{\hat{z}}$ axis. We consider cubic symmetry, that is a finite symmetry group, which in principle is not sufficient for such an angular momentum conservation argument. However, to second order in phonons and magnons, this argument holds anyway. This is because we neglect the effect of $(n^{x})^{2}$, $(n^{y})^{2}$ and $(n^{z})^{2}$, which would produce third order terms in magnons and phonons. That is, $B^{\parallel}$ is not included in the equations, and we are in principle free to set this to any value, as it does not affect the physics. Set it to $B^{\parallel}=B^{\perp}$, and we have in fact assumed isotropic symmetry..

II.0.2 Linearly polarized magnons in biaxial NiO

We will now consider an AFM with two hard-axis anisotropies. For concreteness we will use NiO as an example. To start with, we will neglect its non-trivial spin structure, which introduces the second term in the magnetoelastic Hamiltonian (1). We will therefore solely focus on the first term in the Hamiltonian (1), where only the conventional strain tensor appears. Following the derivation of the Hamiltonian in Appendix C, we find that this is equivalent to assuming that next nearest neighbor interaction is the dominant term contributing to the magnetoelastic interaction.

We start once again from the Hamiltonian (1), which we want to express in terms of the magnon eigenmodes. First, we need to rotate the coordinate system so that the new $\bm{\hat{z}}^{\prime}$ axis coincides with the spin condensation axis. The spins in NiO condense along one of the 12 equivalent $[\bar{1}\bar{1}2]$ directions, within internally ferromagnetic [111] planes. We therefore define the primed coordinate system as $\bm{\hat{z}}^{\prime}=\frac{1}{\sqrt{6}}[-1,\,-1,\,2]$, $\bm{\hat{x}}^{\prime}=\frac{1}{\sqrt{3}}[1,\,1,\,1]$ and $\bm{\hat{y}}^{\prime}=\frac{1}{\sqrt{2}}[-1,\,1,\,0]$. Define the rotation matrix $\mathbf{O}$ so that the primed and un-primed coordinate systems are related by $\bm{r}=\mathbf{O}\bm{r}^{\prime}$. By writing $\bm{n}=\mathbf{O}\bm{n}^{\prime}$, we find

[TABLE]

where we once again have assumed that the magnon and phonons propagate along the $\bm{\hat{z}}$ direction.

In NiO, the magnon eigenmodes are linearly polarized and spinless. In other words, the magnon eigenmodes correspond semi-clasically to oscillations of $n^{\prime x}$ and $n^{\prime y}$ separately. Hence, the Hamiltonian (3) is in fact already given in terms of both the magnon and phonon eigenmodes, and can thus be directly interpreted. We consider first the transverse phonons. Both magnon modes hybridize with both transverse phonon modes, however with different interaction parameters. As a consequence of this, angular momentum in the $\bm{\hat{z}}$ direction is no longer conserved. This is a direct consequence of the lack of rotational symmetry about the $\bm{\hat{z}}$ axis due to the spin condensation axis $\bm{\hat{z}}^{\prime}$ not being aligned with the momentum direction $\bm{\hat{z}}$. Further, we note that only the $n^{\prime x}$ magnon mode hybridizes with the longitudinal phonon mode $\partial u^{z}/\partial r^{z}$. The $n^{\prime x}$ mode is the mode oscillating along the axis with the largest anisotropy, and is followingly the most energetic mode. We therefore conclude that the lower magnon mode passes the longitudinal phonon modes undisturbed, while the upper one is expected to hybridize.

Note that the exact decoupling of the lower magnon mode from the longitudinal phonon modes in the Hamiltonian (3) is a consequence of the magnon eigenmodes being linearly polarized. If the magnon eigenmodes were not linearly polarized, the longitudinal phonon mode would in general couple to both phonon modes. This can be realized by applying an external magnetic field along the $\bm{\hat{z}}^{\prime}$ axis. The effect is that the magnon eigenmodes can be described semi-classically as elliptical precessions of the Néel field around the ground state. For concreteness, let us assume that the eigenmodes are elliptically polarized, $\alpha=\left(An^{\prime x}+iBn^{\prime y}\right)$ and $\beta=\left(Bn^{\prime x}-iAn^{\prime y}\right)$ 222The qualitative results do not depend on the exact form of the eigenmodes as long as they mix $n^{\prime x}$ and $n^{\prime y}$., where $A$ and $B$ depend on the magnetic field strength. The Hamiltonian (3) expressed in terms of the eigenmodes then follows as

[TABLE]

As expected, both magnons $\alpha$ and $\beta$ now hybridize with the longitudinal phonon $\partial u^{z}/\partial r^{z}$.

II.0.3 Antiferromagnets with internally ferromagnetic planes

In the last section, we considered magnon-phonon hybridization in NiO under the assumption that the second term in the Hamiltonian (1) is negligible. We will now look at the effect of the second, spin structure dependent term. If we include only those terms contributing to the magnon-phonon hybridization with momentum along the $\bm{\hat{z}}$ direction 333Additional strain tensor elements appear in the spin structure dependent tensor if we do not require momentum conservation. Additional elements also appear if we consider uniform magnetostriction., the spin structure dependent tensor in NiO is

[TABLE]

The further process of interpreting the hybridization is just the same as shown above; rewrite the Hamiltonian in terms of the magnon eigenmodes, and then read off which modes hybridize. The result is

[TABLE]

Evidently, the spin structure dependent term in the Hamiltonian does not introduce any new types of hybridization, as the magnon and phonon modes which couple are identical to those appearing in Eq. (3). Therefore, the discussion of a simplified NiO-like material in Sec. II.0.2 appears to be valid for the real NiO as well. We predict that the decoupling between a magnon mode and the longitudinal modes can be lifted by applying an external magnetic field along the $\bm{\hat{z}}^{\prime}$ axis. The magnetic field can also be used to smoothly tune the hybridization between the modes, as the coefficients $A$ and $B$ appearing in Eq. (4) depend on the field strength.

III Quantized Hamiltonian

In this section, we will derive the quantized Hamiltonian which is later used to find the exact magnon-polarons in NiO. We start by deriving the magnon Hamiltonian, followed by the phonon Hamiltonian. Last, we will derive the terms which couple magnons and phonons into an effective hybridized state. We stress that although the exact derivations which follow are specific to NiO, the method is fully general and hence valid for all collinear AFMs.

Above the Néel temperature, NiO forms the FCC-structure, whereas it is slightly distorted into a rhombohedral one below Rooksby (1948). This distortion from cubic symmetry is very small, corresponding to an angle of about $0.07^{\circ}$ Slack (1960), and we will therefore neglect it in the following derivations.

III.1 Magnons

The antiferromagnetic ordering in NiO is well established. Below the Néel temperature, the spins order in internally ferromagnetic $(111)$ planes Shull et al. (1951); Hutchings and Samuelsen (1972). The spins on two such neighboring planes are antiparallel, causing the overall structure to be antiferromagnetic. Due to the cubic symmetry of the FCC structure, there are four equivalent $(111)$ planes. In each plane there is moreover a 3-fold degeneracy in the spin direction. There are thus 12 possible antiferromagnetic ground states, one of which are depicted in Fig. 1. As magnons are small spin fluctuations relative to a ground state, we may choose to work from any one of these twelve possible ground states without loss of generalization. We here choose the $(111)$ plane as the plane of internally ferromagnetic order, and $\bm{\hat{z}}^{\prime}=\frac{1}{\sqrt{6}}[-1,\,-1,\,2]$ as the spin quantization axis along which the spins in the (classical) antiferromagnetic ground state are aligned. We moreover let $\bm{\hat{x}}^{\prime}=\frac{1}{\sqrt{3}}[1,\,1,\,1]$ and $\bm{\hat{y}}^{\prime}=\frac{1}{\sqrt{2}}[-1,\,1,\,0]$ define the rest of the primed coordinate system. Note that $\bm{\hat{x}}^{\prime}$ is perpendicular to the $(111)$ plane, and $\bm{\hat{y}}^{\prime}$ is parallel to it. The spins form a bipartite lattice, where the sublattice with spin up (down) is named sublattice A (B).

In the magnetic Hamiltonian we will include exchange interaction and two hard-axis anisotropies. In order to obtain a sufficiently accurate dispersion relation for magnons in NiO, we need to include exchange coupling between both nearest neighbour (nn) and next-nearest neighbour (nnn) spins. Every spin site has 6 nn’s on the same sublattice and 6 nn’s on the opposite sublattice, as well as 6 nnn’s on the opposite sublattice. The magnetic Hamiltonian follows as Lines and Jones (1965); Hutchings and Samuelsen (1972)

[TABLE]

where $\bm{S}_{i}^{\alpha^{\prime}}$ refer to the spin component in the $\alpha^{\prime}$ direction at lattice site $i$, and $D_{x^{\prime}}>0$ and $D_{y^{\prime}}>0$ are anisotropy constants. The summation index $i$ runs over the whole lattice, and $\bm{\delta}_{1}$ and $\bm{\delta}_{2}$ run over nn’s and nnn’s to lattice site $i$, respectively. $J_{1}$ and $J_{2}$ are the corresponding exchange coupling constants. We split the $i$ summations into sums over sublattices A and B with spins $\bm{S}^{\mathrm{A}}$ and $\bm{S}^{\mathrm{B}}$, and do a Holstein-Primakoff transformation of the spin operators in terms of boson operators $a$ and $b$:

[TABLE]

where we have assumed that ${\langle a_{i}^{\dagger}a_{i}^{{\mathstrut}}\rangle/2S\ll 1}$ and ${\langle b_{j}^{\dagger}b_{j}^{{\mathstrut}}\rangle/2S\ll 1}$, and $S=1$ in NiO. We then perform a Fourier transformation of the operators

[TABLE]

where $\bm{x}_{i}$ and $\bm{x}_{j}$ are position vectors on sublattice A and B. Now let $\sum_{\bm{\delta}_{n}\in ab}$ denote the sum over the $n$’th nearest neighbors on sublattice $b$ of a spin belonging to sublattice $a$, where $a,b\in\{\mathrm{A},\,\mathrm{B}\}$. Let $z_{n}^{ab}$ be the number of such neighbors. Use this definition to define the quantity

[TABLE]

One can then show that the Hamiltonian takes the form

[TABLE]

where we have introduced the following coefficients

[TABLE]

This boson Hamiltonian can be diagonalized following the procedure of Ref. [50]. We define

[TABLE]

so that the Hamiltonian takes the form

[TABLE]

Now define $\phi_{\bm{k}}=\mathbf{T}_{\bm{k}}^{-1}\xi_{\bm{k}}$ as the vector of operators that by definition diagonalizes the Hamiltonian. Note that $\mathbf{T}_{\bm{k}}$ in general is a non-unitary transformation, that is $\mathbf{T}_{\bm{k}}^{\dagger}\neq\mathbf{T}_{\bm{k}}^{-1}$. Now define a matrix $\mathbf{g}$ as a commutator between the vector of (bosonic) operators and its Hermitian adjoint, $\mathbf{g}\equiv\commutator{\xi_{\bm{k}}^{{\mathstrut}}}{\xi_{\bm{k}}^{\dagger}}$. By inserting $\xi_{\bm{k}}=\mathbf{T}\phi_{\bm{k}}$ into the commutator, one can show that the tranformation matrix $\mathbf{T}$ satisfies $\mathbf{T}^{\dagger}=\mathbf{g}^{-1}\mathbf{T}^{-1}\mathbf{g}$. The Hamiltonian can then be written into the eigenvalue equation

[TABLE]

where $\mathbf{T}_{\bm{k}i}$ is the $i$’th row of $\mathbf{T}_{\bm{k}}$, and $\hbar\omega_{\bm{k}}$ is the energy of the magnon mode $\bm{k}$. Solving this eigenvalue equation results in

[TABLE]

which describes the dispersion relations for the two free magnon modes in NiO. Experimentally fitted values for the parameters can be looked up in for instance Ref. [48].

III.2 Phonons

NiO forms the FCC structure with two atoms in the basis, one nickel atom and one oxygen atom. These are separated by a distance $a/2$, where $a\approx 4.17\;\AA$ Bartel and Morosin (1971) is the cubic lattice constant. Let $s$ be an index referring to the atom type, with Ni as type 1 and O as type 2. Moreover, let $i$ refer to a site on the FCC lattice at which one such pair of atoms belong. The indices $(i,s)$ therefore uniquely refer to one specific atom in the structure. Now let $u_{is}^{\alpha}$ be the displacement from equilibrium of atom $(i,s)$ in the $\alpha$ direction, where $\alpha\in\{x,y,z\}$. Finally, let $K_{is,jt}^{\alpha\beta}$ be the proportionality constant between the force acting on atom $(i,s)$ in the $\alpha$-direction and the displacement of atom $(j,t)$ in the $\beta$-direction. The phonon energies are then given as the solution of the classical equation of motion

[TABLE]

This equation has $3n$ momentum dependent solutions: $3$ acoustic and $3(n-1)$ optical modes, where $n$ is the number of atoms in the basis. With $n=2$ in NiO, this results in a total of 6 phonon modes. We name these different modes $\omega_{\bm{q}\lambda}$, where $\bm{q}$ is the phonon momentum and $\lambda$ labels the mode. The quantized phonon Hamiltonian follows as

[TABLE]

where $c_{\bm{q}\lambda}^{{\mathstrut}}$ and $c_{\bm{q}\lambda}^{\dagger}$ are phonon annihilation and creation operators.

Eq. (19) is a very general equation, where the phonon energies can be found straightforwardly given a set of force coefficients $K_{is,jt}^{\alpha\beta}$. We will here use the rigid-ion model to compute the phonon eigenmodes in NiO. The rigid-ion model Kellermann (1940); Cochran (1971) is perhaps the simplest model which is able to reproduce a relatively accurate picture of phonons in NiO. It is not the most accurate available model, but it has the advantage of having only a few adjustable parameters, and will suffice for our usage. It is based on modelling each atom as a rigid sphere which moves around an equilibrium position, and is well explained in Refs. [52; 53]. Each atom is connected to its nearby atoms by springs which represent the short-ranged forces between nearby atoms. We include energy terms linear and quadratic in the Ni-O displacements (with force constants $B_{12}$ and $A_{12}$ respectively) and O-O displacements (with force constants $B_{22}$ and $A_{22}$ respectively), while we may neglect the explicit Ni-Ni displacement terms Reichardt et al. (1975).

Additionally, we need to include Coulomb interactions in order to get a realistic model of the optical phonons 444The energy difference between the transverse and longitudinal optical modes depends solely on the effective charge $Z\absolutevalue{e}$. As this energy difference is about 5 THz (in units of $2\pi\hbar$) in NiO, we obviously need to include Coulomb forces in order to attain a realistic phonon dispersion.. Each atom is given an effective charge $\pm Z\absolutevalue{e}$. The Coulomb interaction is long-ranged, meaning that interactions between atoms infinitely far apart contribute. An infinite sum obviously causes numerical difficulties, and to solve this we use a so-called Ewald summation; we split the real space Coulomb summation into a real space sum and a Fourier space integral. We sum over the closest atoms in real space, and approximate the sum over more distant atoms by an integral in Fourier space. This enables us to approximate the formally infinite sum by summing over about $10$ lattice sites in real space. We eventually fit the five constants $A_{12}$, $A_{22}$, $B_{12}$, $B_{22}$ and $Z$ to an experimentally measured phonon dispersion in Ref. [54]. The formal details of this calculation apart from those given above will not be covered here, as it is rather tedious, and we refer the reader to Ref. [52] for further reading.

Each phonon mode is characterized by the atoms moving in a unique pattern. The polarization vector $\bm{\epsilon}_{{\bm{q}\lambda}}$ is formally a six-component vector describing the axes along which the two atom types move for the different phonon modes $\lambda$, as well as the relative phase between the atoms types. As there are three acoustic modes and three optical ones, the only difference between the first and latter three modes is a relative phase $\pi$ between the Ni and O atoms in the optical mode. In the next section, we will use the polarization vector to couple atom displacement and spins. Since only the Ni atoms make a significant contribution to the magnetic Hamiltonian, we will only couple Ni sites. In the following, we may therefore define a three-component polarization vector $\bm{\epsilon}_{{\bm{q}\lambda}}$ describing the movement of the Ni atoms only. As the polarizations of Ni atoms are identical in the acoustic and optical modes, we need only define three different polarizations.

The polarization vector must satisfy the orthogonality relation $\bm{\epsilon}_{\bm{q}\lambda}^{*}\cdotp\bm{\epsilon}_{\bm{q}\lambda^{\prime}}^{{\mathstrut}}=\delta_{\lambda\lambda^{\prime}}$, as well as the completeness relation $\sum_{\bm{q}}\bm{\epsilon}_{\bm{q}\lambda}^{{\mathstrut}}\bm{\epsilon}_{\bm{q}\lambda^{\prime}}^{\dagger}=\mathbb{I}$. Last but not the least, the phonon polarization vectors must be eigenvectors of the equation of motion, Eq. (19). We may conveniently choose the polarization vectors such that $\bm{\epsilon}_{\bm{q}\lambda}^{*}=\bm{\epsilon}_{-\bm{q}\lambda}^{{\mathstrut}}$ Rückriegel et al. (2014). We choose the polarization vectors to be Flebus et al. (2017)

[TABLE]

where $\phi_{\bm{q}}$ and $\theta_{\bm{q}}$ are standard spherical coordinates defining the direction of the momentum $\bm{q}$. $\lambda=1$ and $\lambda=2$ describe transversal modes, while $\lambda=3$ describes longitudinal modes. Recall that these are polarization vectors both for the acoustic and optical modes.

III.3 Magnon-phonon coupling

We will consider magnetoelastic coupling which hybridizes the magnon and phonon modes. Finding the new hybridized eigenstates requires us to diagonalize the Hamiltonian containing magnetic and elastic degrees of freedom. The Hamiltonian under consideration must therefore be quadratic in magnon and phonon operators, meaning that we are only to include interaction terms containing one operator of each sort. The displacement of an ion from equilibrium $\bm{u}_{i}$ is a measure for the elastic degree of freedom, while the spin $\bm{S}_{i}$ at site $i$ is a measure for the magnetic degree of freedom. $\bm{u}_{i}$ is linear in phonon operators (see Eq. (25)), and we may therefore immediately conclude that the interaction term must be linear in $\bm{u}_{i}$.

In Appendix B we do a phenomenological expansion in spins $\bm{S}_{i}$ and displacements $\bm{u}_{j}$ to arrive at two magnetoelastic Hamiltonians. Both terms are linear in displacements $\bm{u}_{j}$, while they are of first and second order in spin. We show that these terms have their origin, among other things, in spin-orbit coupling between a spin and its neighboring ions’ orbital momenta, and in distance dependent exchange interaction, respectively. We will henceforth assume that the term arising from the exchange interaction dominates, and we therefore neglect all terms which are not second order in spin.

The magnetoelastic Hamiltonian under consideration is therefore Evenson and Liu (1969)

[TABLE]

where $i$ is summed over all magnetic lattice sites, $\bm{\delta}$ is a vector pointing from lattice site $i$ to one of its neighboring magnetic atoms, and $\alpha,\beta,\gamma,\lambda\in\{x,y,z\}$ refer to spatial directions. $B_{ij}^{\alpha\beta\gamma\lambda}$ is a tensor of coupling coefficients. $R_{ij}^{\gamma\lambda}$ describes local strains, and we name it the discrete strain tensor. It is defined as

[TABLE]

The discrete strain tensor simplifies to a constant times the continuous strain tensor $\epsilon^{\gamma\lambda}=\frac{1}{2}\left(\frac{\partial u^{\gamma}}{\partial r^{\lambda}}+\frac{\partial u^{\lambda}}{\partial r^{\gamma}}\right)$ in the long-wavelength limit. Note that $R_{ij}^{\alpha\beta}$ is symmetric under exchange of spatial coordinates. We could in principle have coupled the spins to an anti-symmetric elastic tensor as well. An anti-symmetric elastic tensor analogous to $R_{ij}^{\alpha\beta}$ describes local rotations. We will however disregard rotations in our analysis, and couple therefore the spins exclusively to $R_{ij}^{\gamma\lambda}$.

The number of coefficients $B_{ij}^{\alpha\beta\gamma\lambda}$ appearing in the Hamiltonian (22) can be reduced considerably by applying Neumann’s principle, stating that the Hamiltonian must be invariant under symmetry operations of the material itself Birss (1964); Katzir (2004). NiO forms the FCC structure above its Néel temperature, and its structure therefore belongs to the cubic symmetry group $O_{h}$. Neumann’s principle states that the Hamiltonian must be invariant under symmetry operations $\mathcal{R}\in O_{h}$, i.e. ${\mathcal{R}^{-1}\mathcal{H}_{\rm ME}\mathcal{R}=\mathcal{H}_{\rm ME}}$. Furthermore, the Hamiltonian must be translationally invariant. By requiring these symmetries to be fulfilled, we find that the Hamiltonian reduces to

[TABLE]

where $B_{|\bm{\delta}|}^{\alpha\beta}=\delta^{\alpha\beta}B_{|\bm{\delta}|}^{\parallel}+(1-\delta^{\alpha\beta})B_{|\bm{\delta}|}^{\perp}$ Flebus et al. (2017), and $\delta^{\alpha\beta}$ is the Kronecker delta. The $|\bm{\delta}|$ index of $B_{|\bm{\delta}|}^{\alpha\beta}$ means that the coefficients coupling any atoms separated by an equilibrium distance $|\bm{\delta}|$ are equal, which is due to translational and rotational invariance. In other words, there are two coefficients appearing in the Hamiltonian for every $n$’th nearest neighbor spins included in the summation over $\bm{\delta}$.

The displacement vector of a nickel atom at lattice site $i$ can be expressed in terms of the phonon operators as

[TABLE]

where $\omega_{\bm{q}\lambda}$ is the angular frequency of the phonon mode $\lambda$, $m$ is the mass of the nickel atom, and $N$ is the number of nickel lattice sites. The strain tensor between two nickel atoms at position $i$ and $i+\bm{\delta}$ therefore follows as

[TABLE]

If we define the coupling tensor

[TABLE]

the Hamiltonian may be written as

[TABLE]

We now want to expand the Hamiltonian (28) to first order in magnon operators, as this will produce terms quadratic in the boson operators. For this expansion to be justified, we need to expand from the classical ground state of the spins. That is, we first have to express the Hamiltonian in terms of the spins in the primed coordinate system, defined in Sec. III.1. We did a similar procedure in Sec. II.0.2, where we defined $\mathbf{O}$ so that a vector in the primed coordinate system $\bm{r}^{\prime}$ was related to the unprimed coordinates as $\bm{r}=\mathbf{O}\bm{r}^{\prime}$. The spins $\bm{S}$ may then be written as $\bm{S}=\mathbf{O}\bm{S}^{\prime}$. If we then define $\tilde{G}_{{\bm{q}\lambda}\bm{\delta}}^{\alpha\beta}=(\mathbf{O}^{\rm T}\mathbf{G}_{{\bm{q}\lambda}\bm{\delta}}\mathbf{O})^{\alpha\beta}$, the Hamiltonian follows as

[TABLE]

We are now ready to do a Holstein-Primakoff transformation of Eq. (29). We first split the sum over $i$ and $\bm{\delta}$ into four sums: one for each permutation of ${i,i+\bm{\delta}\in\{\mathrm{A},\mathrm{B}\}}$, where A and B are the two sublattices. We will use the following notation: ${\sum_{\bm{\delta}\in ab}}$ means that $\bm{\delta}$ is summed over vectors pointing from a site on a sublattice $a$ to all sites on a sublattice $b$. We note that $\sum_{\bm{\delta}\in\mathrm{AA}}=\sum_{\bm{\delta}\in\mathrm{BB}}$ and $\sum_{\bm{\delta}\in\mathrm{AB}}=\sum_{\bm{\delta}\in\mathrm{BA}}$ due to the equivalency of the sublattices. If we neglect terms which are of third order or higher in the boson operators, and drop the linear terms, the resulting Hamiltonian is

[TABLE]

If we now define

[TABLE]

we may write the Hamiltonian as

[TABLE]

where we have used that $\left(\tilde{G}_{{\bm{q}\lambda}\bm{\delta}}^{\alpha\beta}\right)^{*}=\tilde{G}_{{-\bm{q}\lambda}\bm{\delta}}^{\alpha\beta}$.

We have now expressed the Hamiltonian on a form where magnon and phonon operators are coupled through a single coupling coefficient $\tilde{M}_{{\bm{q}\lambda}}$. All physical details of the material under consideration is contained in $\tilde{M}_{{\bm{q}\lambda}}$. In its definition in Eq. (31), we summed over all neighbors on both sublattices. When performing a calculation, one naturally has to cut off this sum at some point. In NiO, the nn’s of a spin site belonging to sublattice A consist of six sites belonging to sublattice A, and six sites belonging to sublattice B. All nnn’s belong to sublattice B. If we generalize this to more distant neighbors, we find that for any given $n$’th layer of nearest neighbors, if $n$ is odd, then half of the neighbors belong to either sublattice. If $n$ is even, then all neighbors belong to a single sublattice. A natural choice in NiO is therefore to include nearest and next nearest neighbors in the sum over $\bm{\delta}$, as we then include one of each sort of neighbor layers. Two independent sets of magnetoelastic coefficients are therefore included, giving four coefficients in total. We define ${\tilde{B}^{\alpha\beta}=\tilde{B}^{\parallel}\delta^{\alpha\beta}+\tilde{B}^{\perp}\left(1-\delta^{\alpha\beta}\right)}$ as the magnetoelastic coefficient in the nearest neighbor interaction, and ${B^{\alpha\beta}=B^{\parallel}\delta^{\alpha\beta}+B^{\perp}\left(1-\delta^{\alpha\beta}\right)}$ as the magnetoelastic coefficient in the next-nearest neighbor interaction.

IV Magnon-phonon hybridization in NiO

We now combine the magnon Hamiltonian (16), the phonon Hamiltonian (20) and the magnon-phonon Hamiltonian (32). The full Hamiltonian can then be expressed as

[TABLE]

where $\Psi_{\bm{k}}=[\psi^{{\mathstrut}}_{\bm{k}},\,\psi^{\dagger}_{-\bm{k}}]^{\rm T}$ is a vector of all operators, where ${\psi_{\bm{k}}=[a^{{\mathstrut}}_{\bm{k}},\,b^{\dagger}_{-\bm{k}},\,c^{{\mathstrut}}_{1,\bm{k}},\,c^{{\mathstrut}}_{2,\bm{k}},\,c^{{\mathstrut}}_{3,\bm{k}},\,c^{{\mathstrut}}_{4,\bm{k}},\,c^{{\mathstrut}}_{5,\bm{k}},\,c^{{\mathstrut}}_{6,\bm{k}}]}$, and $\mathbf{H}_{\bm{k}}$ is a $(16\times 16)$ non-diagonal matrix. We diagonalize $\mathbf{H}_{\bm{k}}$ with the same procedure as we did with the magnon Hamiltonian in Section III.1. This diagonalization reveals the new hybridized eigenstates with the corresponding energy eigenvalues.

In order to get an initial overview of the full momentum dependence of the energy dispersions, we have plotted the energies of the free magnons and free phonons in the first phonon Brillouin zone in Fig. 2. The black lines depict the phonons, while the red lines depict the magnons. Note that the magnons are non-degenerate, which is due to the hard-axis anisotropies. As can be observed in the figure, the magnon modes cross the optical phonon modes in two distinct areas, at frequencies about $11.3$ THz and $17.3$ THz. We will refer to these areas as the first and second crossing point, respectively. These are the areas where the magnon-phonon hybridization becomes apparent, and where the modes are neither magnon- nor phonon-like. Where the hybridization is strong, the properties mix, and the modes should rather be labeled magnon-polarons. Note that the modes do not cross near the zone center, as they do in for instance YIG Kikkawa et al. (2016); Flebus et al. (2017).

We will now include magnetoelastic coupling in the analysis. We want to display the results for realistic values of the magnetoelastic coefficients. However, to the best of our knowledge, the magnetoelastic coefficients in NiO are not precisely determined. There exist magnetostriction measurements Slack (1960) which in principle may be used to determine magnetoelastic coefficients, but neither of these measurements are sufficiently detailed to determine all four coefficients we use in this analysis. We will therefore rather assess the expected order of magnitude of the coefficients. We do this by assuming that all but one coefficient are negligible, and use the magnetostriction measurements presented in Ref. [46] to estimate the remaining coefficient. As a result, we find that the coefficients take values between approximately [math] THz and $100$ THz. Given this approximative method, we do not expect the following results to be quantitatively accurate apart from the order of magnitude. However, we expect the qualitative effect of each coefficient to be accurate. Combined with the descriptive Eqs. (3) and (6), we are able to supply a thorough analysis of the qualitative magnon-phonon hybridizations in NiO. This may in turn easily be generalized to other cubic collinear AFMs.

We have plotted the magnon-phonon dispersion in the (anti-)crossing areas in Figs. 3 and 4 for different values of the magnetoelastic coefficients. As the modes now mix, the previous black/red labeling for phonons/magnons can no longer be applied. In this and all following figures, all modes are thus colored differently in order for them to be easily recognized. We have continued the assumption from the approximate assessment of the coefficients, namely that we assume that all but one coefficient are negligible, and therefore display the (anti-)crossings with only one non-zero coefficient at a time. We display the modes for three different non-zero values of the magnetoelastic coefficients in the range which was found to be realistic: 0, 50 and 100 THz.

Fig. 3 displays the (anti-)crossings between the magnons and the TO phonon modes. There are a number of features in this plot which should be addressed. First, $B^{\parallel}$ causes no hybridization between the modes. This is simply due to that $B^{\parallel}$ only couples to longitudinal phonon modes, and thus does not affect the TO phonons. Second, $\tilde{B}^{\parallel}$ and $B^{\perp}$ apparently cause all modes to hybridize. The same conclusion can be drawn by directly reading off the hybridizations from Eqs. (3) and (6). Third, the lower magnon mode does not couple to any TO phonons for any values of $\tilde{B}^{\perp}$. Looking at Eq. (6), we find that the magnon mode associated with fluctuations in $n^{\prime y}$ does not couple to the TO phonon modes if only $\tilde{B}^{\perp}$ is non-zero. As $\bm{\hat{y}}^{\prime}$ is the axis with the weakest hard axis anisotropy, $n^{\prime y}$ corresponds to the lowest energy magnon, and this therefore confirms the result of Fig. 3.

The (anti-)crossings between the magnons and the LO phonon mode are displayed in Fig. 4. All qualitative features of this plot may be explained by analyzing Eqs. (3) and (6). First, $\tilde{B}^{\parallel}$ and $B^{\perp}$ cause no hybridization of the modes. This follows directly from the semi-classical equations, as these coefficients do not couple to $\partial u^{z}/\partial r^{z}$. Second, $\tilde{B}^{\perp}$ and $B^{\parallel}$ makes the upper magnon mode hybridize with the LO phonon, also in line with the predictions of Eqs. (3) and (6). Third, the lower magnon mode do not couple to the longitudinal phonon mode at all. Having recognized the lower magnon mode as oscillations in $n^{\prime y}$, this result was also implied by the semi-classical analysis.

Another prediction of the semi-classical analysis in Secs. II.0.2 and II.0.3 was that turning on a magnetic field would lift the decoupling of the LO phonon mode and the lower magnon mode. We therefore supply additional plots of the (anti-)crossings between the magnon modes and the LO phonon mode, this time with an external magnetic field present along the $\bm{\hat{z}}^{\prime}$ axis, in Fig. 5. We do this for three different magnetic field strengths, 0, 1 and 2 T, well below the spin-flop field where the quantum theory is expected to be imprecise to this order in the magnon operators Machado et al. (2017). The four plots display the results for the four permutations of $\tilde{B}^{\parallel},B^{\perp}\in\{0,\,25\}$ THz. The plots confirm that applying a magnetic field couples all modes. Increasing the magnetic field strength shows hybridization between the previosuly uncoupled modes, which confirms that the hybridization is tunable. We have limited the magnetoelastic coefficients to 25 THz simply because this gives more readable plots. Larger $B^{\parallel}$ causes the hybridization to be very strong, as displayed in the lower left plot of Fig. 4. This causes the upper magnon mode to cross the lower magnon mode at lower $k$, and the effect of applying a magnetic field would hence be most evident at lower $k$, making the plots somewhat less coherent. This effect can be seen in the two lower plots, where $B^{\parallel}$ is non-zero. The qualitative effect of applying a magnetic field is nonetheless also present for larger values of the magnetoelastic coefficients.

V Discussion

The quantitative analysis in the previous section has provided results for the magnon-phonon hybridizations in NiO given a selection of values for the magnetoelastic coefficients. We have furthermore found that the semi-classical analysis given in Secs. II.0.2 and II.0.3 is able to describe all of these hybridizations qualitatively. Both approaches are fairly general, and the methods can straightforwardly be applied to other cubic collinear AFMs. The specific results for NiO can also to some extent be generalized to other cubic collinear antiferromagnetic systems, at least qualitatively.

Antiferromagnetic cubic materials introduce at least two new features to magnon-polarons as compared to their ferromagnetic counterparts. The first feature is that longitudinal phonon modes hybridize with magnons in many antiferromagnetic structures. There are two origins of this hybridization: first, the longitudinal modes hybridize if the anisotropies cause the spins to condense non-parallel to any of the crystal axes; second, this occurs due to the second term in the magnetoelastic Hamiltonian (1), which is there due to non-trivial spin ordering in the ground state. The first origin is not unique to AFMs. Indeed, the spins condense in the [111] directions in ferromagnets such as pure nickel and magnetite Bozorth and Walker (1950); Medrano et al. (1999). One should expect the longitudinal phonons propagating along the crystal axes to hybridize with magnons in these ferromagnetic materials as well. This follows directly from a ferromagnetic analogy of Eq. (4), where the Néel field is substituted with magnetization and where one of the magnon polarizations is discarded (the evanescent mode). The second origin however, is only attainable if there are at least two spin sublattices involved, as this is the least requirement for creating a non-trivial spin-structure. This latter effect introduce hybridization between longitudinal phonons and magnons even if the spins are aligned along one of the crystal axes in the ground state.

The second feature we have observed is that the hybridization between antiferromagnetic magnons and the phonons is tunable by an external magnetic field. Moreover, the hybridization may be switched off entirely in certain structures, and thus shows potential for binary control. This is a unique feature of bipartite magnetic structures. Although magnon-phonon hybridizations in ferromagnets can be tuned in the sense that the magnon energy can be increased relative to the phonon modes, the interaction in itself cannot be completely turned on and off. This may however be done in certain AFMs, where NiO is a prominent example of a material in which this phenomenon should occur. This is due to there being two magnon modes in AFMs, which enables "moving" the magnetoelastic interaction between the magnon modes by changing the properties of the eigenstates. More precisely, since applying a magnetic field changes the spin oscillations of the magnon eigenmodes qualitatively, and because the hybridizations is sensitive to this oscillation pattern, we may use the magnetic field to tune the hybridization in AFMs.

We expect magnon-polarons in AFMs to be found at optical phonon energies, which are usually well up in the THz range. In contrast, magnon-polarons in ferromagnets are typically found at the crossings between magnons and acoustic phonons. This is due to the linear dispersion of antiferromagnetic magnons at low $k$. The effect of this is that the magnon dispersions never cross the acoustic phonon dispersions if their velocity is greater than that of the phonons. Followingly, there is no strong hybridization between the magnons and acoustic phonons. The high frequencies at which magnons-polarons are found in AFMs make the accessibility of antiferromagnetic magnon-polarons lower than the ferromagnetic ones. This in turn might reduce their applicability. For instance, magnon-polarons have shown to play an important role in the SSE effect in YIG Kikkawa et al. (2016); Flebus et al. (2017), as the properties of phonons affect the induced spin current. This would not occur in most AFMs, as the magnon-polarons are not necessarily thermally accessible below the Néel temperature of the material.

An important point to address is, precisely, how to access antiferromagnetic magnon-polarons. Both crossings between magnon-like and phonon-like modes in NiO occur far away from the zone center. This stands in contrast to ferromagnetic magnon-polarons, which are typically found at low $k$. For that reason, magnon-polarons are not accessible by for instance conventional first order Raman scattering, which aims to excite modes at very long wavelengths due to the negligible momentum of the photon. We have already discarded thermal excitation as an alternative, due to the high energies. The magnon-polarons may be accessed with neutron scattering, which previously has been used to map the dispersion relations of both magnons and phonons separately in NiO. The most promising way of accessing the magnon-polarons in AFMs might however be with femtosecond optics Kampfrath et al. (2010); Ju et al. (2004); Battiato et al. (2010); Berritta et al. (2016).

Injection of coherent phonons at high frequencies in the THz range has been achieved using ultrafast lasers Ruello and Gusev (2015). In a direct analogy with spin pumping driven via coherent elastic waves Weiler et al. (2012), and subsequently formed magnon-polarons Kamra et al. (2015); Hayashi and Ando (2018), in ferromagnets, these coherent THz phonons may directly excite the corresponding magnon-polarons, and result in a spin pumping current Kamra et al. (2015); Hayashi and Ando (2018). The latter may be detected electrically via inverse spin Hall effect (ISHE). Since a magnetic field may be used to tune the hybridizations, one might even be able to tune the measured ISHE voltage with the magnetic field.

VI Conclusion

We explore magnon-polarons in cubic collinear antiferromagnets, and focus on their qualitative difference to their ferromagnetic counterparts. We find that there are two distinct features of antiferromagnetic magnon-polarons. First, antiferromagnetic materials with either a complex spin structure or spin alignment axis non-parallel to any of the cubic axes generally result in hybridization between magnons and longitudinal phonons. Second, the hybridizations may be tuned by an external magnetic field by changing the qualitative properties of the magnons eigenmodes. NiO is an example of such a material, where a hybridization may even be turned on and off.

Magnons-polarons in antiferromagnets seem to have a reduced applicability compared to ferromagnetic ones due to their high energies and locations at intermediate momenta. However, the rapid evolving fields of ultrafast dynamics and femtosecond optics provide tools for an easier access to the antiferromagnetic magnon-polarons. We suggest spin pumping driven via coherently injected phonons as a promising approach for investigating these magnon-polarons.

Acknowledgements

This work was supported by the Research Council of Norway through its Centres of Excellence funding scheme, Project No. 262633 "QuSpin", as well as by the European Research Council via Advanced Grant No. 669442 "Insulatronics".

Appendix A Derivation of the magnetoelastic Hamiltonian

Magnetoelastic coupling relates magnetic and elastic degrees of freedom. In a discrete lattice, at a site $i$, the magnetic moment is proportional to the spin $\bm{S}_{i}$. A measure of the elastic degree of freedom is the displacement of atom $i$, $\bm{u}_{i}$. The magnetoelastic coupling depends on the spin $\bm{S}_{i}$ and displacement $\bm{u}_{j}$ at all sites $i$ and $j$. We will now discuss phenomenological models of the magnetoelastic coupling starting with the simplest possible forms. We consider materials in which the magnetic atoms form a Bravais lattice, all of which are invariant under inversion symmetry.

Let us begin by discussing a simple ansatz, that the discrete magnetoelastic model is bilinear in spin and displacement,

[TABLE]

Here, $B_{ij}^{\alpha\beta}$ is a phenomenological coupling tensor relating the spin at site $i$ with the atom displacement at site $j$. Neumann’s principle states that the physical properties of a crystal must share the symmetries of the crystal Birss (1964); Katzir (2004). Under a transformation $\mathcal{R}$, the spin transforms as ${\bm{S}\rightarrow\absolutevalue{\mathcal{R}}\mathcal{R}\bm{S}}$ and the displacement transforms as ${\bm{u}\rightarrow\mathcal{R}\bm{u}}$. Since we have assumed that the magnetic atoms form a Bravais lattice, the Hamiltonian is invariant under the inversion operation $\mathcal{R}=-1$ and $\absolutevalue{\mathcal{R}}=-1$. Consequently, $B_{ij}^{\alpha\beta}=-B_{ij}^{\alpha\beta}$. The only solution is $B_{ij}^{\alpha\beta}=0$. Hence, the Hamiltonian (34) does not contribute.

Let us proceed by first restoring invariance under the inversion operation without expanding to higher orders in the magnetic or elastic degrees of freedom. One may observe that introducing another quantity that transforms as a vector fulfills the requirement of invariance. In our system, the only natural vector we have left is the position vector $\bm{r}_{i}$. Our next attempt is therefore

[TABLE]

which will have finite elements $B_{ijk}^{\alpha\beta\gamma}$ even when requiring inversion invariance.

There are also other constraints. The Hamiltonian (35) must be invariant under uniform translations of the lattice. Mathematically, this can be expressed by a uniform displacement $u_{j}^{\beta}\rightarrow u_{j}^{\beta}+\delta u^{\beta}$, or by a uniform shift of the position vectors $r_{k}^{\gamma}\rightarrow r_{k}^{\gamma}+\delta r^{\gamma}$. The Hamiltonian (35) must be invariant under both of these transformations separately.

We begin by considering a uniform displacement, where the invariance requirement reads

[TABLE]

This relation must hold for all different spin configurations $\{S_{i}^{\alpha}\}$ and atom configurations $\{r_{k}^{\gamma}\}$. Furthermore, as $\delta u^{\beta}$ is arbitrary, the relation must hold for every component $\beta$. The resulting constraint for the magnetoelastic coefficients is therefore

[TABLE]

We now make use of the second invariance requirement. A uniform shift of the position vectors $r_{k}^{\gamma}\rightarrow r_{k}^{\gamma}+\delta r^{\gamma}$ should leave the Hamiltonian unchanged:

[TABLE]

This relation must hold for any spin configuration $\{S_{i}^{\alpha}\}$ and atom displacements $\{u_{j}^{\beta}\}$. As the shift $\delta r^{\gamma}$ is arbitrary, the relation must hold for every component $\gamma$. Thus, we are left with the constraint

[TABLE]

Requiring the Hamiltonian to be translationally invariant induces two constraints on the tensor $B_{ijk}^{\alpha\beta\gamma}$, Eqs. (37) and (39).

Let us now inspect the Hamiltonian (35) more closely. $S_{i}^{\alpha}$ and $u_{j}^{\beta}$ are dynamical variables. In contrast, $r_{k}^{\gamma}$ is fixed when the lattice properties are defined. In other words, $\{r_{k}^{\gamma}\}$ defines the equilibrium lattice, and the dynamics related to displacements from equilibrium is contained in $\{u_{j}^{\beta}\}$. Keeping this in mind, we could define an effective coupling tensor as ${\tilde{B}_{ij}^{\alpha\beta}=\sum_{k}\sum_{\gamma}B_{ijk}^{\alpha\beta\gamma}r_{k}^{\gamma}}$, so that the Hamiltonian (35) reads

[TABLE]

At first sight Eq. (40) might seem to have become analogous to the starting ansatz (34), the latter of which does not contribute. However, there is an important distinction since $\tilde{B}_{ij}^{\alpha\beta}$ is not a tensor of constant coefficients. Instead, $\tilde{B}_{ij}^{\alpha\beta}$ is a sum of products between a tensor of constant coefficients and position vector components. Therefore, the transformation properties of $\tilde{B}_{ij}^{\alpha\beta}$ differ from the ones of $B_{ij}^{\alpha\beta}$. As a result, the Hamiltonian (40) essentially differs from the starting Hamiltonian (34).

By expressing the Hamiltonian (35) as in (40), we realize that several $k$ and $\gamma$ components of $B_{ijk}^{\alpha\beta\gamma}$ contribute to the effective tensor $\tilde{B}_{ij}^{\alpha\beta}$, but the relative contribution to the sum is not important. The only physical significance of the introduction of the position $r_{k}^{\gamma}$ is its transformation properties. This implies that, without a loss of generality, we can choose a selection of the tensor elements $B_{ijk}^{\alpha\beta\gamma}$ to be equal to zero as long as we do not break any symmetries of the lattice under consideration. There is an infinite number of such choices in an infinite lattice. We will follow a path that is physically transparent because, in the continum limit, it couples the spins to strain tensor components. As will be evident, we obtain this by leaving $B_{ijk}^{\alpha\beta\gamma}$ finite for the following indices $i$, $j$, and $k$: if $i\neq j$ then $k\in\{i,j\}$, and if $i=j$ then $k$ can point to any lattice site. $B_{ijk}^{\alpha\beta\gamma}$ is set to zero for all other $k$’s.

We may now insert the definition of the non-zero tensor elements $B_{ijk}^{\alpha\beta\gamma}$ into the constraints (37) and (39). The constraint (37) then reads

[TABLE]

The second constraint (39) reads

[TABLE]

In all of these equations, $i\neq j$. Eqs. (41) and (43) imply ${B_{iij}^{\alpha\beta\gamma}=B_{iji}^{\alpha\beta\gamma}=-B_{ijj}^{\alpha\beta\gamma}}$. We therefore define a new tensor ${B_{ij}^{\alpha\beta\gamma}\equiv B_{ijj}^{\alpha\beta\gamma}}$, and insert all the constraints above into the Hamiltonian (35). The resulting Hamiltonian is

[TABLE]

where both $i$ and $j$ run over all lattice sites. This Hamiltonian (45) is thus the lowest order non-zero magnetoelastic Hamiltonian.

The Hamiltonian (45) describes both rotations, which are antisymmetric in $\beta\gamma$, and strains, which are symmetric in $\beta\gamma$. We can separate these two effects by defining the matrices

[TABLE]

where $R^{\beta\gamma}_{\mathrm{S},ij}$ captures strains and $R^{\beta\gamma}_{\mathrm{R},ij}$ captures rotations. In the following, we will restrict the analysis to strains only, and we therefore drop the rotation term. We name $R^{\beta\gamma}_{\mathrm{S},ij}$ the discrete strain tensor, and denote it simply as $R^{\beta\gamma}_{ij}$ from now on. By summing over the vector $\bm{\delta}\equiv\bm{r}_{j}-\bm{r}_{i}$ instead of $j$, we may rewrite the Hamiltonian as

[TABLE]

The $ij$-index of $B_{ij}^{\alpha\beta\gamma}$ was changed to $|\bm{\delta}|$. This is possible due to translational and rotational invariance of the Bravais lattice. Due to the normalization factor $|\bm{r}_{i}-\bm{r}_{j}|^{-2}$ in the definition of the strain tensor (46), the coupling tensor elements appearing in (48) are related to the tensor elements in (45) by the inverse of this normalization factor. As a result, $B_{0}^{\alpha\beta\gamma}=0$, and the first non-zero coefficients appear in the nearest neighbor interaction.

We may now use Neumann’s principle to derive the selection rules of $B_{|\bm{\delta}|}^{\alpha\beta\gamma}$. We find that there are 18 independent coefficients in triclinic crystals, 8 in monoclinic crystals, and fewer as we increase the symmetry. In cubic crystals, we find that $B_{|\bm{\delta}|}^{\alpha\beta\gamma}=0$ is the only solution. We conclude that the Eq. (48) cannot describe any magnetoelastic coupling in cubic crystals.

As we want the magnetoelastic Hamiltonian to be able to describe cubic antiferromagnets, we must include an additional term. Our next attempt is quadratic in the spin degrees of freedom,

[TABLE]

Just as we did for $\mathcal{H}^{\rm II}_{\rm ME}$, we require the Hamiltonian to be invariant under uniform translations of the lattice, that is ${u_{k}^{\gamma}\rightarrow u_{k}^{\gamma}+\delta u^{\gamma}}$ and ${r_{l}^{\lambda}\rightarrow r_{l}^{\lambda}+\delta r^{\lambda}}$, and obtain the constraints

[TABLE]

In the following, we will use a local approximation. We assume that the interaction between $\bm{S}_{i}$, $\bm{S}_{j}$ and $\bm{u}_{k}$ is dominated by the terms where $k\in\{i,j\}$. Intuitively, this follows if we view the Hamiltonian (49) as a distance dependent exchange interaction. Hence, the relevant displacements are the displacements of the involved spins. In other words, we can view it as a local expansion in the lattice distortions around the spins.

As above, where we discovered an arbitrariness in the indices of the coupling tensor related to the position vector $r_{l}^{\lambda}$, we may choose to set $B_{ijkl}^{\alpha\beta\gamma\lambda}=0$ for a selection of $l$’s without a loss of generality. We would like the Hamiltonian (49) to be consistent with the continuum limit result of spins coupling to the strain tensor. We therefore choose $B_{ijkl}^{\alpha\beta\gamma\lambda}$ to be non-zero only if $l\in\{i,j\}$, in addition to the already mentioned $k\in\{i,j\}$. We now insert the choice of non-zero tensor elements into the constraints of Eqs. (50) and (51) derived above. The first constraint (50) then reads

[TABLE]

Additionally, the second constraint (51) reads

[TABLE]

In all constraints, $i\neq j$. The constraints (52) and (55) together imply $B_{ijii}^{\alpha\beta\gamma\lambda}=-B_{ijij}^{\alpha\beta\gamma\lambda}=-B_{ijji}^{\alpha\beta\gamma\lambda}$, and the constraints (53) and (56) additionally imply $B_{ijii}^{\alpha\beta\gamma\lambda}=B_{ijjj}^{\alpha\beta\gamma\lambda}$. Let us therefore define $B_{ij}^{\alpha\beta\gamma\lambda}\equiv B_{ijii}^{\alpha\beta\gamma\lambda}$. The resulting Hamiltonian follows as

[TABLE]

If we consider strains only, and disregard rotations, the Hamiltonian reads

[TABLE]

where we changed the summation variable $j$ to $\bm{\delta}=\bm{r}_{j}-\bm{r}_{i}$ just as we did in the derivation of (48). $\bm{\delta}$ runs over all lattice sites except $\bm{\delta}=0$, fulfilling constraint (54).

One may verify that the Hamiltonian (58) has non-zero contributions both for cubic crystals and for crystals subject to a uniform strain. We conclude that the Hamiltonians in Eqs. (48) and (58) combined give the lowest order phenomenological picture of magnetoelastic coupling in magnetic crystals forming a Bravais lattice.

Appendix B Physical origin of the magnetoelastic Hamiltonians

The purpose of this section is to discuss the physical origins of the magnetoelastic Hamiltonians (48) and (58) that we derived in Appendix A. To elucidate the properties, we use as a starting point well-known interactions, and show how these generate Eqs. (48) and (58). First, we consider spin-orbit interaction between a spin and the orbit of its neighboring ions and demonstrate that this leads to a Hamiltonian equivalent to (48). Second, we consider a distance dependent exchange interaction and observe that this reproduces the Hamiltonian (58).

Consider first the spin-orbit coupling between a spin at site $i$ and the orbital magnetic momentum of an ion at site $j$. The magnetic moment of the spin is $\bm{\mu}^{\rm spin}_{i}=\gamma_{\rm s}\bm{S}_{i}$, where $\gamma_{\rm s}$ is the gyromagnetic ratio of the spin. The orbital magnetic moment of ion $j$ in the rest frame of spin $\bm{S}_{i}$ is $\bm{\mu}^{\rm ion}_{j}=\gamma_{\rm ion}\bm{L}_{j}$, where $\gamma_{\rm ion}$ is the gyromagnetic ratio of the ion, and $\bm{L}_{j}$ is its orbital angular momentum. The orbital angular momentum is ${\bm{L}_{j}=m(\bm{r}_{i}-\bm{r}_{j})\crossproduct\partial_{t}(\bm{u}_{i}-\bm{u}_{j})}$, where $m$ is the mass of the ion and $\partial_{t}=\partial/(\partial t)$ is the time differential operator. A general form of the corresponding spin-orbit Hamiltonian follows as

[TABLE]

where all constants are contained in the coupling tensor $C_{ij}^{\alpha\beta\gamma}$. Assuming plane wave solutions of the displacements, ${u_{i}^{\gamma}=U_{\bm{k}}^{\gamma}\exp{i\bm{k}\cdot\bm{r}_{i}-i\omega t}}$, where $\bm{k}$ is the wave vector of the plane wave, gives $\partial_{t}(u_{i}^{\gamma}-u_{j}^{\gamma})=-i\omega(u_{i}^{\gamma}-u_{j}^{\gamma})$. The resulting spin-orbit Hamiltonian becomes

[TABLE]

where $\tilde{C}_{ij}^{\alpha\beta\gamma}(\omega)=-i\omega C_{ij}^{\alpha\beta\gamma}$. We recognize that the expression (60) is equivalent to the Hamiltonian (48) after symmetrization in $\beta\gamma$.

We next consider a general exchange interaction (not excluding Dzyaloshinskii-Moriya interaction),

[TABLE]

where $A_{ij}^{\alpha\beta}$ is a tensor of coefficients coupling spin components $\alpha$ and $\beta$ at lattice sites $i$ and $j$. $A_{ij}^{\alpha\beta}$ depends on the distance separating the spins $\bm{S}_{i}$ and $\bm{S}_{j}$. Let $\bm{r}_{i}-\bm{r}_{j}$ be the equilibrium position vector separating the spins, and let $\bm{u}_{i}-\bm{u}_{j}$ be the displacement relative to equilibrium. We may then expand $A_{ij}^{\alpha\beta}$ to first order in the relative displacement as

[TABLE]

If we insert the expansion (62) into the Hamiltonian (61), we see that a new term coupling the spins at site $i$ and $j$ to the relative displacement appears. This interaction is proportional to the coupling tensor $\partial A_{ij}^{\alpha\beta}/\partial(u_{i}^{\gamma}-u_{j}^{\gamma})\rvert_{\bm{u}_{i}-\bm{u}_{j}=0}$. Now note that this coupling tensor is not a tensor of constant coefficients, like $A_{ij}^{\alpha\beta}$. It transforms differently due to the operator $\partial/\partial(u_{i}^{\gamma}-u_{j}^{\gamma})$ working on it. For instance, under inversion, ${\partial/\partial(u_{i}^{\gamma}-u_{j}^{\gamma})\rightarrow-\partial/\partial(u_{i}^{\gamma}-u_{j}^{\gamma})}$. For convenience, we therefore introduce the coupling tensor $B_{ijkl}^{\alpha\beta\gamma\lambda}$, implicitly defined by

[TABLE]

where $k$ and $l$ run over all lattice sites. The components of $B_{ijkl}^{\alpha\beta\gamma\lambda}$ transform just as scalars under inversion, and we therefore choose to proceed with $B_{ijkl}^{\alpha\beta\gamma\lambda}$ as the coupling tensor. As before, we note that $(r_{k}^{\gamma}-r_{l}^{\gamma})$ is not a dynamical variable, but merely a constant with the desired transformation properties once the lattice has been defined. Its introduction does therefore not alter the physical content of the Hamiltonian.

We find the non-zero tensor elements $B_{ijkl}^{\alpha\beta\gamma\lambda}$ by requiring the Hamiltonian (61) to respect the symmetries of the lattice. Following the argument presented in Appendix A, we may pick only the $B_{ijkl}^{\alpha\beta\gamma\lambda}$’s where $k,l\in\{i,j\}$ to be non-zero. We define $B_{ij}^{\alpha\beta\gamma\lambda}\equiv B_{ijii}^{\alpha\beta\gamma\lambda}$, and may then write the exchange Hamiltonian to first order in the relative displacement as

[TABLE]

We recognize that the second term is equivalent to the Hamiltonian (58) after symmetrization in the $\gamma\lambda$ indices in order to disregard rotations.

Appendix C Long-wavelength magnetoelastic Hamiltonian

We will now derive the long-wavelength magnetoelastic Hamiltonian for a cubic collinear antiferromagnets, starting from the general form given in Eq. (24). In a collinear antiferromagnet, there are two sublattices, A and B. In the (classical) ground state, the spins on each sublattice are anti-parallel. Motivated by the bipartite lattice, we first separate the sums over $i$ and $\bm{\delta}$ in Eq. (24) into two separate contributions, one in which $i$ is a site on sublattice A, and one in which $i$ is a site on sublattice B. We can further separate each of these sums into two sums, where the vector $\bm{\delta}$ points between sites on the same sublattice, or between the sublattices. In total, we thus have four separate sums. However, because sublattice A and B are equivalent, there are only two independent sums. We therefore have

[TABLE]

where $\bm{\delta}\in$ AA (AB) denotes that $\bm{\delta}$ points from a site on sublattice A to a site on sublattice A (B). We will now transit to the long-wavelength limit. First define two spin fields, $\bm{S}^{\rm A}(\bm{r})$ and $\bm{S}^{\rm B}(\bm{r})$, living on sublattice A and B respectively. We then introduce the Néel vector $\bm{n}(\bm{r})=\frac{1}{2}\left(\bm{S}^{\rm A}(\bm{r})-\bm{S}^{\rm B}(\bm{r})\right)$ and the local magnetization $\bm{m}(\bm{r})=\frac{1}{2}\left(\bm{S}^{\rm A}(\bm{r})+\bm{S}^{\rm B}(\bm{r})\right)$. We will in the following neglect the local magnetization, and may therefore express the Hamiltonian as

[TABLE]

Notice the relative minus sign between the sum over AA and AB, which has been introduced because the spins on the two sublattices are approximately anti-parallel, $\bm{n}=\bm{S}^{\rm A}(\bm{r})\approx-\bm{S}^{\rm B}(\bm{r})$. This sign change is characteristic to antiferromagnets and strongly affects the final result in certain materials.

In order to fully transit to the long-wavelength limit, we first need to do the sum over $\bm{\delta}$. In order to do so, we expand the Néel vector to first order,

[TABLE]

For notational simplicity, we define $\left.\frac{\partial n^{\alpha}(\bm{r})}{\partial r^{\gamma}}\right\rvert_{\bm{r}=\bm{r}_{i}}\equiv D^{\alpha\gamma}(\bm{r_{i}})$. The Hamiltonians then reads

[TABLE]

Now rearrange the sums over $\bm{\delta}$ into separate sums arising from the different layers of $n$’th nearest neighbors at distances $|\bm{\delta}_{n}|$. That is, write

[TABLE]

Note that we have put $B_{\absolutevalue{\bm{\delta}}}^{\alpha\beta}$ outside the sums, as these coefficients are equal for all neighbors included in an $n$’th nearest neighbor summation. We now need to evaluate only two kinds of sums, namely $\sum_{\bm{\delta}_{n}}R_{i,i+\bm{\delta}_{n}}^{\alpha\beta}$ and $\sum_{\bm{\delta}_{n}}R_{i,i+\bm{\delta}_{n}}^{\alpha\beta}\delta_{n}^{\gamma}$.

In cubic collinear antiferromagnets, there are three types of $n$’th nearest neighbor layers to any given spin $i$ on sublattice A: 1. All neighbors belong to sublattice A and are therefore parallel to $\bm{S}_{i}$. 2. All neighbors belong to sublattice B, and are therefore anti-parallel to $\bm{S}_{i}$. 3. Half of the spins belong to sublattice A and the other half belong to sublattice B. In situation 1 and 2, all $n$’th nearest neighbors belong to one sublattice. We will now assume that only the $B_{\absolutevalue{\bm{\delta}_{n}}}^{\alpha\beta}$’s where $|\bm{\delta}_{n}|$ is much shorter than the wavelength of lattice strains contributes to the sum. In other words, we assume that the range of the magnetoelastic interaction is relatively short. In cubic crystals, we then have to good approximation $\sum_{\bm{\delta}_{n}}R_{i,i+\bm{\delta}_{n}}^{\alpha\beta}=C_{n}\epsilon^{\alpha\beta}(\bm{r}_{i})$ where $C_{n}$ is some constant, and $\sum_{\bm{\delta}_{n}}R_{i,i+\bm{\delta}_{n}}^{\alpha\beta}\bm{\delta}^{\gamma}_{n}=0$. Situation 1 and 2 thus couple elements of the Néel field straightforwardly to elements of the conventional strain tensor.

In situation 3, however, half of the terms in the summation over $\bm{\delta}_{n}$ comes with a minus sign, as is evident from Eq. (68). By evaluating the sums over $\bm{\delta}$ arising from situation 3 explicitly, we find $\sum_{\bm{\delta}_{n}}^{\prime}R_{i,i+\bm{\delta}_{n}}^{\alpha\beta}=\tilde{C}_{n}\tilde{\epsilon}^{\alpha\beta}(\bm{r}_{i})$ where the primed summation indicates that terms arising from a $\bm{\delta}_{n}$ pointing between the sublattices are accompanied with a minus sign, $\tilde{C}_{n}$ is some constant, and $\tilde{\epsilon}^{\alpha\beta}(\bm{r}_{i})$ is some tensor with elements that are linear combinations of strain tensor elements. In general, the structure of this tensor depends on the microscopic spin structure of the material. The last sum gives no contribution for situation 3 either, $\sum_{\bm{\delta}_{n}}^{\prime}R_{i,i+\bm{\delta}_{n}}^{\alpha\beta}\bm{\delta}^{\gamma}_{n}=0$, meaning that the gradient of the Néel field does not contribute to the final Hamiltonian to this order. Situation 3 hence couples the Néel field to another tensor $\tilde{\epsilon}$, not equal to the strain tensor. This tensor only arises when there exist nearest neighbor layers in which a portion of the spins are part of sublattice A and a portion are part of sublattice B. NiO is an example of such a material, and likewise are CoO, FeO and MnO.

We may now finalize the transition into the long-wavelength limit, by taking ${\sum_{i}\rightarrow\int\mathrm{d}\bm{r}}$. The magnetoelastic Hamiltonian for cubic collinear antiferromagnets then follows as

[TABLE]

where $B^{\alpha\beta}=\sum_{\mathrm{s}(n)\in\{1,2\}}C_{n}B_{|\bm{\delta}_{n}|}^{\alpha\beta}$ and $\tilde{B}^{\alpha\beta}=\sum_{\mathrm{s}(n)\in 3}\tilde{C}_{n}B_{|\bm{\delta}_{n}|}^{\alpha\beta}$, where the notation $\mathrm{s}(n)\in\{1,2\}$ means summing over the $n$’th nearest neighbor layers belonging to either situation 1 or 2 as introduced in the above paragraph, and $\mathrm{s}(n)\in 3$ means summing over the $n$’th nearest neighbor layers belonging to situation 3. Each tensor contains two independent coefficients, $B^{\alpha\beta}=B^{\parallel}\delta^{\alpha\beta}+B^{\perp}\left(1-\delta^{\alpha\beta}\right)$ and $\tilde{B}^{\alpha\beta}=\tilde{B}^{\parallel}\delta^{\alpha\beta}+\tilde{B}^{\perp}\left(1-\delta^{\alpha\beta}\right)$.