The full magnon spectrum of yttrium iron garnet

The magnetic insulator yttrium iron garnet can be grown with exceptional quality, has a ferrimagnetic transition temperature of nearly 600 K, and is used in microwave and spintronic devices that can operate at room temperature. The most accurate prior measurements of the magnon spectrum date back nearly 40 years, but cover only 3 of the lowest energy modes out of 20 distinct magnon branches. Here we have used time-of-flight inelastic neutron scattering to measure the full magnon spectrum throughout the Brillouin zone. We find that the existing models of the excitation spectrum fail to describe the optical magnon modes. Using a very general spin Hamiltonian, we show that the magnetic interactions are both longer-ranged and more complex than was previously understood. The results provide the basis for accurate microscopic models of the finite temperature magnetic properties of yttrium iron garnet, necessary for next-generation electronic devices. Complete spin-wave dispersions have been measured on a crystal of Yttrium Iron Garnet (YIG) by inelastic neutron scattering. Andrew Princep and co-workers from Oxford University (in partnership with Innovent e.V) have grown high quality YIG crystals and measure the dispersion of magnons, which are wave-like quantized collective excitations of spins, in up to energies of 120 meV using inelastic neutron scattering. They further reproduce the entire spectrum by using linear spin–wave theory, considering various forms of the magnetic exchange interactions and pathways. The analysis confirms the important role of nearest-neighbour magnetic exchange interactions, but calls for reinterpreting the nature of longer-ranged interactions. These results offer essential information on the optical magnon modes at room temperature and are of importance for technological applications of YIG.


INTRODUCTION
Yttrium iron garnet (YIG) is the 'miracle material' of microwave magnetics. 1 Since its synthesis by Geller and Gilleo in 1957, 2 it is widely acknowledged to have contributed more to the understanding of electronic spin-wave and magnon dynamics than any other substance. 3 Magnons are bosonic quasiparticles and are the quanta of magnetic oscillations of systems of periodically ordered magnetic moments. YIG (chemical formula Y 3 Fe 5 O 12 , crystal structure depicted in Fig. 1a) is a ferrimagnetic insulating oxide with the lowest magnon damping of any known material. Its exceptionally narrow magnetic resonance linewidth-orders of magnitude lower than the best polycrystalline metals-allows magnon propagation to be observed over centimetre distances. This makes YIG both a superior model system for the experimental study of fundamental aspects of microwave magnetic dynamics 4 (and indeed, general wave and quasi-particle dynamics 5,6 ), and an ideal platform for the development of microwave magnetic technologies, which have already resulted in the creation of the magnon transistor and magnon logic gates. 4,7,8 The unique properties of YIG have underpinned the recent emergence of new fields of research, including magnonics: 9 the study of magnon dynamics in magnetic thin-films and nanostructures, 4,10 and magnon spintronics, concerning structures and devices that involve the interconversion between electronic spin currents and magnon currents. 1 Such systems exploit the established toolbox of electron-based spintronics as well as the ability of magnons to be decoupled from their environment and efficiently manipulated both magnetically and electrically. 8,[11][12][13] There has been particular interest in YIG after the demonstrated transmission of electrically injected magnons over distances of 1 mm. 5 Short-wavelength magnons can also diffuse over distances of~50 μm even at room temperature, ideal for application in devices. 14 Significant interest has also developed in quantum aspects of magnon dynamics, resulting in the Bose-Einstein condensation of magnons even at room temperature, 6 or using YIG as the basis for new solid-state quantum measurement and information processing technologies including cavity-based QED, optomagnonics, and optomechanics. 15 It has also recently been realised that one can stimulate strong coupling between the magnon modes of YIG and a superconducting qubit, potentially as a tool for quantum information technologies. 16 Spin caloritronics has also recently emerged as a potential application of YIG, utilising the spin Seebeck effect (SSE) and the spin Peltier effect (SPE) to interconvert between magnon and thermal currents, either for efficient large-scale energy harvesting, or the generation of spin currents using thermal gradients. 17 If the research into classical and quantum aspects of spin wave propagation in YIG is to achieve its potential, it is absolutely clear that the community requires the deep understanding of its mode structure which only neutron scattering measurements can offer. Since magnons are bosonic, the calculation of their thermal population and transport properties (such as the magnon thermal conductivity and the SSE) at non-zero temperatures is both intrinsically quantum-mechanical and dependent on the precise distribution of magnon modes.
In many theories and experiments, YIG is treated as a ferromagnet with a single, parabolic spin wave mode, 18,19 simply because the influence of YIG's complex electronic and magnetic structure on spin transport is not known in sufficient detail. Such approaches must break down at temperatures when the optical modes are appreciably populated, and a detailed knowledge of the structure of the optical modes is a necessary first step in any realistic model of the magnetic properties of YIG in this operational regime. Despite this, surprisingly little data exists relating to the detail of its magnon mode structure. The previous key work in this area is due to Plant et al., 20,21 35-45 years ago. Using a triple-axis spectrometer, these early neutron scattering measurements were able to record a few of the spin wave modes up to approximately 55 meV (14 THz), but crucially there are a total of 20 distinct modes for a given wave vector and they are predicted to extend up to approximately 90 meV (22 THz). 3 There is a large cluster of poorly understood modes in the region of 30-50 meV which are known to soften substantially (of order 30%) by room temperature, owing largely to the reduction in sublattice magnetisation. This places many of these modes at energies where they are expected to have an appreciable thermal population at 300 K, and so their precise distribution must be known to enable calculations of thermal magnetic properties.

RESULTS
Data were collected (see Methods) as a large, 4-dimensional hypervolume in energy and momentum space, covering the complete magnon dispersion over a large number of Brillouin zones. Figure 2 shows two-dimensional energy-momentum slices from this hypervolume with the wave vector along three highsymmetry directions, with scattered intensities normalised to a measurement on vanadium (see Methods). A large number of modes can be seen up to an energy of 80 meV, whilst data in other slices show modes extending up to nearly 100 meV. The spectrum is dominated by a strongly dispersing and well-isolated acoustic mode at low energies (the so-called 'ferromagnetic' mode), and a strongly dispersing optical mode separated from this by a gap of approximately 30 meV at the zone centre. Intersecting this upper mode is a large number of more weakly dispersing optical modes in the region of 30-50 meV.
The magnetic order in YIG is driven by exchange interactions, which are the result of the quantum mechanical requirement that the electronic states are antisymmetric under particle exchange. We model these using a Heisenberg effective spin Hamiltonian, which is appropriate as YIG is both a good insulator and the Fe 3+ ions (S = 5 / 2 , g = 2) possess a negligible magnetic anisotropy due to the quenched orbital moment: where the summation is over pairs of spins. As a check, we also evaluated models which included anisotropic terms in H but found these terms to be vanishingly small, consistent with previous results. The spin Hamiltonian was diagonalized using the SpinW software package 22 and the calculated magnon dispersion was fitted by a constrained nonlinear least squares method to 1D cuts taken through the 2D intensity slices (see Supplementary  Information). We do not include any scaling factors for the magnon intensity, so the agreement between the model and the data in terms of absolute intensity is indicative of the quality of the model. Our final/best-fit model includes isotropic exchange interactions up to the 6th nearest neighbour, labelled J 1 -J 6 in Fig. 1b. The parameters in this work can be mapped to the exchanges commonly considered for YIG as follows: J ad = J 1 , J dd = J 2 , and J aa = {J 3a , J 3b }, where the subscript refers to the majority tetrahedral (d) and minority octahedral (a) sites. The inclusion of interactions up to the 6th nearest neighbour is guided by previous experimental results and ab initio calculations which indicate that the magnitudes of J 4 -J 6 could be as much as 5-10% that of the dominant exchange interaction J 1 . 21,23 Due to the extremely large number (20) of magnetic atoms within the primitive cell, and the consideration of so many exchange pathways, this analysis would be impossible without the use of sophisticated software such as SpinW as the construction of an analytic model would be prohibitively time consuming. During the fitting process, features in the spectrum were weighted so that weak but meaningful features in the data were considered as significant as strong features. The SpinW model output is then convoluted with the calculated experimental resolution of the MAPS spectrometer, including all features of the neutron flight path and associated focussing/defocussing effects, as well as the detector coverage and effects from symmetrisation (see Supplementary Information for details). The final fitted values of the exchanges are listed in Table 1, alongside values obtained from the two studies by Plant, and the calculations performed in recent ab initio work.

DISCUSSION
An important difference between our results and those of previous authors 23 is that there are two symmetry-distinct 3rdnearest-neighbour bonds (the so-called J aa in the literature, which we label J 3a and J 3b ) which have identical length but differ in symmetry. The J 3b exchange lies precisely along the bodydiagonal of the crystal, and thus has severely limited symmetryallowed components owing to the high symmetry of the bond  . The J 3a exchange connects the same atoms with the same radial separation, but represents a different Fe-O-Fe exchange pathway as a result of the different point group symmetry (point group C 2 ), so it is distinguished from J 3b by the environment around the Fe atoms. Models including anisotropic exchange or Dzyaloshinskii-Moriya interactions on the 1st-4th neighbour bonds were tested, but such interactions were found to destabilise the magnetic structure for arbitrarily small perturbations. We also find that J 2 (J dd ) is much smaller than previously supposed-the main effect of this exchange is to increase the bandwidth and split the optic modes clustered around 40 meV in a way contradicted by the data.
It has been pointed out 24 that the magnetic structure of YIG is incompatible with the cubic crystal symmetry, although to date no measurements have found any evidence for departures from the ideal cubic structure. Nevertheless, it is necessary to refine the magnetic structure in a trigonal space group (a symmetry that is experimentally observed in terbium rare earth garnets where the magnetoelastic coupling is much stronger 25 ), in order to obtain a satisfactory goodness of fit, and a magnetic moment which agrees with bulk magnetometry. 26 Treating the unit cell of YIG in this fashion for the purposes of the SpinW simulation would be feasible, but would introduce a large number of free parameters which (given the very small size of the departure from cubic symmetry) would nevertheless be expected to change very little from a cubic model. This expectation is borne out by the excellent agreement between the data and a spin wave model that assumes cubic symmetry (see Fig. 2

and the Supplementary Information for more details).
Features absent from the data which are relevant to technological applications (such as the conversion of microwave photons into magnons) include any strong indications of magnon-magnon or magnon-phonon coupling. The temperature of the experiment (10 K) results in a very low magnon thermal population, strongly suppressing the 2-and 4-magnon processes which are necessary for magnon-magnon effects to be apparent. However, the model we have produced would allow a full, quantum mechanical, multi-magnon calculation to be done which could reproduce multi-magnon effects at higher temperatures. A strong magnonphonon coupling would be expected to cause both broadening and anomalies in the dispersion of the magnon modes. 27 We do not observe any such effects, although our measurements would not be sensitive to any magnon-phonon coupling that shifts or broadens the spin wave signal by less than about 3 meV (the instrumental resolution). This absence of detectable magnon-phonon coupling is not unexpected given the extremely long spin-lattice relaxation times observed in YIG. 6 Our results require a substantial revision of the impact of the optical modes on the room-temperature magnetic properties. These differences are compared with the existing work in Fig. 3, where we plot the antisymmetric combination of transverse spectral functions S xy (Q,ω) − S yx (Q,ω) averaged over many Brillouin zones, which is proportional to the sign and magnitude of the measured spin-Seebeck current arising from the associated magnon modes. 28 Results are plotted for both models in the centred (i.e. non-primitive) cell. Compared with the previous model used to calculate the spin-Seebeck current there is a The exchange interactions are defined in Fig. 1, and are per pair of spins Fig. 3 Simulations of the spectrum dispersion in YIG. a The previous model 3, 28 and b our new model of the magnon dispersion, where the colour and intensity correspond to the sign and magnitude of the spectral function S xy (Q,ω) − S yx (Q,ω) averaged over many Brillouin zones, which is responsible for the spin Seebeck effect. 28 Results are plotted for the centred (i.e. non-primitive) cell. Red and blue modes have positive and negative values of the spectral function and are associated with spin currents that flow in opposite directions. The horizontal dashed line indicates k B T at room temperature compression and shift of the 'positively' polarised (red) optical modes (i.e. those modes which would precess counterclockwise with respect to an applied field) in the 30-50 meV region as well as a redistribution of spectral weight. As has recently been shown, 28 the thermal population and softening of the 'negative' (blue) modes at elevated temperatures substantially reduces the magnitude of the measured spin-Seebeck current. A full calculation of the temperature dependent magnon population in a multimagnon picture, based on our model, can be used to derive strict limitations on device performance and determines the optimum operating temperature as well as an exact calculation of the magnon thermal conductivity which is very difficult to access experimentally. 29,30 Our results show that the distribution of optical modes is very different from what had previously been assumed, which also has consequences for the temperaturedependent softening and broadening of the modes. 28 Our new measurements can therefore be used as the basis for a precise microscopic model of the temperature dependent dynamical magnetic properties in YIG.
We have also estimated the parabolicity of the lowest lying 'ferromagnetic' magnon mode, i.e. the point where the deviation from a quadratic fit becomes greater than 5%. We find this region to extend~10% of the way towards the Brillouin zone boundary along the Γ-H direction (the (0,0,1) direction in the centred unit cell), representing about 0.1% of the entire Brillouin zone. A quadratic fit to the model dispersion yields a spin-wave stiffness constant D = ħω/k 2 of 532 meV Å 2 (85.2 × 10 -41 Jm 2 ), in good agreement with previous measurements. 31 The departure from a purely parabolic acoustic magnon dispersion, as well as the population of optic magnon modes directly generates the temperature dependence of the spin Seebeck effect and our model can be used to fully understand such effects even at elevated temperatures through extension to a multi-magnon picture following the procedure in refs. 28, 32. We have presented the most detailed and complete measurement of the magnon spectrum of YIG in a pristine, high quality crystal. Our work has uncovered substantial discrepancies between previous models and the measured dispersion of the optical magnon modes in the 30-50 meV range, as well as the total magnon bandwidth and the detailed nature of the magnetic exchanges. Using linear spin-wave theory analysis and through a detailed consideration of possible forms of the exchange interactions and the symmetries of the exchange pathways, we are able to reproduce the entire magnon spectrum across a large number of Brillouin zones including a quantitatively accurate reproduction of the absolute neutron scattering intensity. We confirm the importance of the nearest-neighbour exchange, but are forced to radically reinterpret the nature and hierarchy of longer-ranged interactions. Our findings are of importance to technological applications of YIG, particularly those utilising the spin Seebeck effect or the magnon thermal conductivity, which are very sensitive to the optical magnons in the region of 30-50 meV. The results overturn 40 years of established work on magnons in YIG, and the detailed spin Hamiltonian we have presented allows any magnetic property of YIG to be calculated as a function of temperature. We foresee this to be particularly useful for modelling and understanding magnon transport phenomena at room temperature.

Crystal growth
YIG crystal growth was carried out in high-temperature solutions applying the slow cooling method. 33 Starting compounds of yttrium oxide (99.999%) and iron oxide (99.8%) as solute and a boron oxide-lead oxide solvent were placed in a platinum crucible and melted in a tubular furnace to obtain a high-temperature solution. 34 Using an appropriate temperature gradient only a few single crystals nucleate spontaneously at the cooler crucible bottom and forced convection, obtained by accelerated crucible rotation technique (ACRT), allows a stable growth which results in nearly defect-free large YIG crystals. 35 The YIG crystal used in this study exhibits a size of 25 mm × 20 mm × 11 mm and a weight of 12 g. It was confirmed by neutron and X-ray diffraction that the YIG crystal was a single grain with a crystalline mosaic of approximately 0.07 degrees FWHM. Preliminary measurements were made on a crystal that was grown by the optical floating-zone method, starting from a pure powder of YIG.
Neutron scattering data collection and reduction Data were collected on the MAPS time-of-flight neutron spectrometer at the ISIS spallation neutron source at the STFC Rutherford Appleton Laboratory, UK. On direct geometry spectrometers such as MAPS, monochromatic pulses of neutrons are selected using a Fermi chopper with a suitably chosen phase. In our experiment neutrons with an incident energy (E i ) of 120 meV were used with the chopper spun at 350 Hz, giving energy resolution of 5.4 meV at the elastic line, 3.8 meV at an energy transfer of 50 meV, and 3.1 meV at an energy transfer of 90 meV. The spectra were normalized to the incoherent scattering from a standard vanadium sample measured with the same incident energy, enabling us to present the data in absolute units of mb sr -1 meV -1 f.u. -1 (where f.u. refers to one formula unit of Y 3 Fe 5 O 12 ). Neutrons are scattered by the sample onto a large area detector on which their time of flight, and hence final energy, and position are recorded. The two spherical polar angles of each detector element, time of flight, and sample orientation allow the scattering function S(Q, ω) to be mapped in a four dimensional space (Q x , Q y , Q z , E). In our experiment the sample was oriented with the (HHL) plane horizontal, while the angle of the (00L) direction with respect to the incident beam direction was varied over a 120°range in 0.25°steps. This resulted in coverage of a large number of Brillouin zones, which was essential in order to disentangle the 20 different magnon modes. Due to the complex structure factor resulting from the number of Fe atoms in the unit cell, the mode intensity varies considerably throughout reciprocal space. The large datasets recording the 4D space of S(Q,ω),~100 GB in this case, were reduced using the Mantid framework, 36 and both visualized and analyzed using the Horace software package. 37 Taking advantage of the cubic symmetry of YIG, the data were folded into a single octant of reciprocal space (H > 0, K > 0, L > 0), adding together data points that are equivalent in order to produce a better signal-to-noise. 2D slices were taken from the 6 reciprocal space directions depicted in Fig. 2 and in Supplementary Information.

Data availability
The neutron scattering data presented in this work were recorded on the MAPS spectrometer at the ISIS Facility, and are available from https://doi. org/10.5286/isis.e.73945114.