Field-angle resolved magnetic excitations as a probe of hidden-order symmetry in CeB6

P. Y. Portnichenko,1 A. Akbari,2, 3, 4, 5, ∗ S. E. Nikitin,1, 2 A. S. Cameron,1 A. V. Dukhnenko,6, † V. B. Filipov,6 N. Yu. Shitsevalova,6 P. Čermák,7, 8 I. Radelytskyi,7 A. Schneidewind,7 J. Ollivier,9 A. Podlesnyak,10 Z. Huesges,11 J. Xu,11 A. Ivanov,9 Y. Sidis,12 S. Petit,12 J.-M. Mignot,12 P. Thalmeier,2 and D. S. Inosov1, ‡ 1Institut für Festkörperund Materialphysik, Technische Universität Dresden, 01069 Dresden, Germany 2Max-Planck-Institut für Chemische Physik fester Stoffe, Nöthnitzer Str. 40, 01187 Dresden, Germany 3Asia Pacific Center for Theoretical Physics, Pohang, Gyeongbuk 790-784, Korea 4Department of Physics, POSTECH, Pohang, Gyeongbuk 790-784, Korea 5Max Planck POSTECH Center for Complex Phase Materials, POSTECH, Pohang 790-784, Korea 6I. M. Frantsevich Institute for Problems of Material Sciences of NAS, 3 Krzhyzhanovsky Street, 03680 Kiev, Ukraine 7Jülich Center for Neutron Science at MLZ, Forschungszentrum Jülich GmbH, Lichtenbergstraße 1, 85748 Garching, Germany 8Faculty of Mathematics and Physics, Department of Condensed Matter Physics, Charles University, Ke Karlovu 5, 121 16, Praha, Czech Republic 9Institut Laue-Langevin, 71 Avenue des Martyrs CS 20156, 38042 Grenoble Cedex 9, France 10Neutron Scattering Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA 11Helmholtz-Zentrum Berlin für Materialien und Energie GmbH, Hahn-Meitner-Platz 1, 14109 Berlin, Germany 12Laboratoire Léon Brillouin, CEA-CNRS, CEA/Saclay, 91191 Gif sur Yvette, France (Dated: January 7, 2020)


I. INTRODUCTION
The high degeneracy and strong Coulomb repulsion of f-electrons in lanthanide and actinide compounds can lead to exotic quantum matter states at low temperatures [1][2][3][4][5][6][7]. The hybridization of valence (conduction) electrons with strongly correlated f-electrons may result in the formation of heavy-fermion metals with quasiparticles that have large effective masses and opening of hybridization gaps that can lead to the Kondo insulator state. Furthermore, at even lower temperatures, broken-symmetry phases (usually magnetic order or superconductivity) may appear driven by residual quasiparticle interactions. Of particular interest are "hidden order" (HO) states which consist in spontaneous order of higher rank (rank r ≤ 2l = 6) multipoles of the f electrons (l = 3). Here even and odd r correspond to the preservation and breaking of time-reversal symmetry by the order parameter, respectively [2]. In this language, charge or valence (monopole) ordering and spin (dipole) ordering correspond to rank 0 and rank 1 order parameters, respectively. The most well-known examples of materials that host higherrank order are cubic CeB 6 (rank 2 quadrupole and rank 3 octupole) [8,9] and tetragonal URu 2 Si 2 (proposed rank 5 dotriakontapole [10][11][12][13] or rank 4 hexadecapole [14]). Further examples are the cubic 4 f-skutterudites [3], 1-2-20 cage compounds [15], and R 3 Pd 20 X 6 (R = rare earth, X = Si, Ge) clathrates [16][17][18][19][20]. The symmetry of HO states is difficult to identify because the common x-ray and neutron diffraction (in zero field) yield no signal. Alternative techniques like neutron diffraction in external field [21], resonant x-ray scattering [22][23][24], or ultrasonic investigations [25,26] had to be applied to identify combined AFQ and field-induced AFO order in phase II of CeB 6 below T Q = 3.2 K.
Additional methods for characterizing HO are desirable. One way is to look at the magnetic excitation spectrum in the HO phase which carries the imprint of the multipolar interactions and the HO parameter in its dispersion relations [6,27]. In the conventional approach, the dispersion is calculated using a specific candidate model for the HO in some fixed applied field and then compared to the dynamical structure factor measured by single-crystal inelastic neutron scattering (INS). This approach was also tried for CeB 6 where theoretical results, employing extended Holstein-Primakoff (HP) [28] and random-phase approximation (RPA) methods [29][30][31] have been used to calculate the dispersion and intensity of the magnetic excitation spectrum. Early [32] and later more accurate [33,34] experiments showed that a few modes at high-symmetry points in the simple-cubic Brillouin zone (BZ) could be identified, however an overall determination of the dispersion in the whole BZ still remains elusive. This is not surprising due to several reasons. First of all, to minimize the number of adjustable parameters in the theoretical models, they have been so far restricted to only nearest-neighbor interactions among the multipoles. However, generalized Ruderman-Kittel-Kasuya-Yosida (RKKY) interactions that are mediated by the conduction electrons are expected to be long-range with an oscillatory character in direct space [35][36][37]. There is plenty of experimental evidence that RKKY-type exchange dominates the multipolar interactions in CeB 6 , as it was shown that the diffuse momentum-space distribution of quasielastic magnetic intensity in the high-temperature paramagnetic state nicely matches the static electronic response (Lindhard) function of the conduction electrons [38,39]. Recent calculations of the effective RKKY-type exchange terms between different types of multipoles in CeB 6 [36,37] demonstrate that secondnearest-neighbor contributions are always stronger than the nearest-neighbor ones, and even third-nearest-neighbor terms are not negligible. According to these results, truncating the interactions at ∼10% of the dominant one would still result in a model with at least 5 independent couplings that would have to be treated as adjustable parameters when fitting the experimental magnon dispersions along the conventional approach. To the best of our knowledge, so far no spin-dynamical calculations in the AFQ state that would conform to the complexity of this realistic multi-parameter model have been attempted. The good news, however, is that individual further-neighbor interactions should mainly influence the dispersion away from the zone center, whereas zone-center magnetic excitations are determined by the total sum of all neighbor interactions for each pair of sublattices. For the latter case, only the total sum enters in the Zeeman energy due to the molecular field and in the RPA denominator of the collective magnetic response in Eq. (6). Therefore zone-center modes should be much less sensitive to the various individual coupling constants between Ce 3+ ions. We therefore expect that magnetic excitations at the Γ point should show better agreement with the currently available models of multipolar excitations. Moreover, zone-center excitations can be probed not only by INS, but also by electron spin resonance (ESR) [40][41][42][43], hence the experimental results at the Γ point can be obtained from two complementary spectroscopic techniques [43] and are therefore more reliable than elsewhere in the BZ.
The second problem in comparing the experimental INS data with the theoretical models consists in the proximity to the antiferromagnetic (AFM) dipolar phase III that represents the ground state of CeB 6 at low magnetic fields, which is not captured by the previous HP and RPA approaches. Moreover, the intensity of magnetic excitations near highsymmetry Γ and R points in the simple-cubic BZ is enhanced by conduction-4 f electron Kondo screening, which leads to itinerant heavy 4 f quasiparticles [44] and ultimately to spin-exciton formation within the AFM phase [33,45]. This possibility is not contained in the theory based on completely localized 4 f-electron picture, hence a reliable comparison of experimental data with the present models is only possible at high magnetic fields deep within phase II, where the AFM order parameter is suppressed.
The third difficulty is that the calculated inelastic magnetic response in phase II is limited to the dipolar response function, whereas strictly speaking neutron scattering is sensitive to all odd-rank magnetic multipoles (such as dipole, octupole, dotriakontapole etc.). At short scattering vectors, |Q| → 0, the form factor is expected to suppress all higherorder multipolar contributions, which justifies the so-called dipolar approximation that has been employed for calculating the dynamical structure factor even for multipolar-ordered systems. The non-monotonic form factor of higher-order multipoles vanishes at Q = 0 and then starts to increase until reaching a maximum at some finite momentum transfer. This has been established from both theory and experiment for elastic neutron scattering [3,[46][47][48], but to the best of our knowledge, the highly involved theory of INS beyond the dipolar approximation [49,50] was never successfully applied to calculate the dynamical response functions for any compound with a multipolar-ordered ground state. On the other hand, experiments clearly demonstrate that the intensity of the Γ -point excitation in CeB 6 increases with wave vector when going from the shortest Γ (001) to the secondshortest Γ (110) reciprocal-lattice vector [34], suggesting a non-monotonic form factor characteristic of multipolar moments that cannot be described in the dipolar approximation. Access to the same excitation at even longer wave vectors is restricted by the worsening of energy resolution imposed by the kinematic constraints, and therefore the actual share of non-dipolar spectral weight in the experimental INS spectrum remains unknown. We expect that it should change with the applied magnetic field, as it induces secondary dipolar and octupolar order parameters on top of the primary AFQ order, thereby activating additional magnetic degrees of freedom that are subject to collective fluctuations. The multipolar corrections should increasingly influence the measured intensities of magnetic excitations towards higher |Q|, leading to deviations from the established spin-dynamical models.
In high magnetic fields, sufficiently far above phase III, when the Kondo and spin-exciton effects due to itineracy of heavy f-electrons are suppressed, the localized approach becomes appropriate. In particular, the gapped Goldstone mode around the Γ point, which was predicted by the theory [30,31], has been eventually observed in an INS experiment [34], in quantitative agreement with the ferromagnetic resonance seen earlier using ESR spectroscopy [40][41][42][43]. However, due to the difficulties outlined above, we found no quantitative agreement between the theory and the observed dispersion of multipolar excitations throughout the BZ. In the present work, we propose an alternative approach to analyze the fingerprints of HO in the magnetic excitation spectrum, which appears to be more promising in providing quantitative information while staying within the same theoretical framework. Rather than changing continuously the momentum transfer and following the dispersion of multipolar excitations in momentum space, one may fix the wave vector at some high-symmetry point with large INS intensity (e.g. R or Γ ) and instead change the direction and strength of the applied magnetic field. Ideally, the field could be continuously rotated, so that changing the field alone would give access to a three-dimensional parameter space for every single value of the momentum-transfer vector Q. However, due to the limited measurement time, we were practically restricted only to several high-symmetry directions of the magnetic field in this work.
We demonstrate that not only the energy and intensity, but even the number of experimentally observable multipolar excitations changes as a function of field direction, offering an abundant source of experimental information that remained virtually unexplored in all the previous measurements. Thermodynamically, the anisotropy of the critical fields of phase II in CeB 6 was found to be determined by the underlying AFQ/AFO hidden order. Therefore, a similar anisotropy in the pattern of collective modes under the field rotation is not at all surprising and can be calculated within the framework of available models. We will demonstrate the usefulness of this new approach to HO symmetry on the canonical example of CeB 6 multipole model, where the order parameter is already known from other methods. Our work has the twofold aim of proposing a general new method for HO investigation with INS for fixed momentum transfer in a rotating field and making concrete experimental and theoretical analysis of magnetic excitation modes in CeB 6 as a function of field direction and magnitude, in particular as guidance for further experiments and future theoretical developments.
The paper is organized as follows. In Sec. II we describe new INS experiments on CeB 6 obtained in fields up to 16.5 T. We present data measured in several high-symmetry points of the BZ for different field directions, obtained in multiple INS experiments on the same single-crystalline sample of CeB 6 . Then, in Secs. III A and III B we recapitulate the basic RPA-type response theory for calculating the magnetic excitation spectrum for a localized Γ 8 model in the presence of quadrupolar and induced octupolar HO. In Sec. III C we present the results for the multipolar mode dispersions using composite Γ + 5 /Γ − 2 HO model of CeB 6 as visible in the intensity of the dipolar structure factor. In Sec. III D we present analogous calculations in dependence on the field strength and direction in the BZ high-symmetry points. We then give a detailed comparison of the theoretical and experimental results, showing clear evidence for a new high-energy branch of multipolar excitations observed for the first time in CeB 6 . Finally, we will summarize our results in Sec. IV, where we argue that our results on CeB 6 serve as a proof of principle for a new alternative way of obtaining quantitative information about hidden-order symmetry and interactions between multipole degrees of freedom in f -electron materials with non-dipolar order parameters. This method should extend the conventional approach to the analysis of neutronscattering data and may find applications in a much broader class of correlated-electron systems.

A. Sample description and experimental configurations
Most of the experimental results presented in this paper were measured on the same rod-shaped Ce 11 B 6 single crystal with a mass of ∼4 g, which was previously used for experiments reported in Refs. [33,34,43]. To minimize neutron absorption by the highly absorbing 10 B isotope present with a 19.9% abundance in natural boron, the crystal was grown by the floating-zone method from the isotope-enriched 11 B powder (Ceradyne Inc., 99.6 % enrichment level), as described elsewhere [33,see Methods]. The growth axis of the crystalline rod was approximately parallel to the [110] crystallographic axis, with a deviation of a few degrees. Before every experiment, the crystal was pre-aligned in the corresponding scattering plane (which varied from one experiment to another) using an x-ray Laue camera.
We have measured the magnetic excitation spectrum of CeB 6 by neutron spectroscopy for different field directions using 6 different instruments. These are the coldneutron triple-axis spectrometers (TAS) PANDA (Heinz  [60,61]. The latter were required to measure the magnetic excitation spectrum at a constant magnetic field over the entire BZ. The corresponding experimental configurations and sample orientations are summarized in Table I. In every TAS experiment, a cold Be filter was used in order to avoid higher-order neutron contamination from the monochromator, and the final wave vector of the neutrons was fixed to k f = 1.3 Å −1 (E f = 3.50 meV) or 1.5 Å −1 (4.66 meV), as indicated for every dataset. TOF measurements were done with the incident neutron wavelength fixed at λ i = 5 Å (E i = 3.27 meV) for IN5 and 5.1 Å (3.15 meV) for CNCS experiments. The magnetic field was applied in  . The data are offset vertically for clarity. They illustrate the appearance of resonance C (smaller peak) at a higher energy with respect to resonance A (stronger peak) at magnetic fields above 14 T. Both peak positions shift upwards in energy with increasing field. Solid lines are Lorentzian fits. the vertical direction (perpendicular to the scattering plane) using various cryomagnets available at the corresponding neutron facilities. For PANDA experiments, we used either the 7.5 T liquid-free vertical-field cryomagnet with a 3 He insert or the 5 T and 8 T LHe-cooled vertical-field cryomagnets for measurements with magnetic field along the [110], [001] and [111] crystal directions, respectively. For the experiment at 4F2 we used a vertical-field 9 T magnet. In the FLEXX and ThALES experiments, we used similar 15 T magnets equipped with a lambda point refrigerator. To achieve the maximal field of 16.5 T, the same magnet was used with an additional dysprosium "booster", which concentrates magnetic field lines in a smaller sample volume, so that higher magnetic field can be achieved at the expense of the smaller sample size. For this unique configuration, a smaller piece of an identical crystal rod (∼ 0.9 g in mass) had to be cut to fit into the sample space of the "booster". In the IN5 and CNCS experiments, we used the low-background 2.5 T "orange"-type cryostat based magnet and the 5 T cryomagnet, respectively.

Zone-center excitations
We start the presentation of our results with the high-field INS data obtained with the magnetic field B [110], because this field direction has been most extensively studied in the past, and in particular in our earlier work [43] we already presented both the field dependence for fields up to 7 T and the momentum dependence of multipolar excitations for certain values of the field in this direction. We found perfect quantitative agreement between the energy of the single resonance at the Γ point and the ferromagnetic resonance observed in ESR (resonance A) within the measured field range. At the same time, a second ESR resonance (resonance B) was seen at fields above 12 T [41], which we did not cover in our initial experiment. We therefore decided to extend our INS measurements up to higher fields, expecting to find the emergence of the resonance B at lower energies relative to the main resonance, as suggested by the ESR data. However, the corresponding energy scans measured with setup 1 at the base temperature of 1.7 K, which are presented in Fig. 1, show no signatures of such an additional peak in the expected energy range of 0.9-1.1 meV up to the maximal field of 14.5 T accessible with this particular experimental setup. The background in our range of interest is clean and essentially constant, and the energy resolution is sufficiently narrow so that the tail of the peak corresponding to resonance A does not overlap with the expected position of resonance B at such high magnetic fields. For instance, at 14.5 T the expected peak splitting is 0.35 meV, while the instrument resolution is less than 0.1 meV [52]. Therefore, we can conclude that resonance B is not visible to neutrons, possibly due to certain selection rules that are different from those of ESR. To confirm this unexpected result, we repeated the experiment at the high-flux spectrometer ThALES (setup 3) [55], the corresponding data are presented in Fig. 2 (a). The absence of the peak associated with resonance B was reproduced, but extending the spectrum to higher energies revealed a new peak at a higher energy of ∼2.25 meV in the 14.5 T dataset that rapidly vanished when the field was decreased below its maximal value. This result suggests that an additional resonance (resonance C) forms above ∼14 T, at approximately the same field as resonance B, but on the opposite side of the intense peak corresponding to resonance A. Because 14.5 T was the highest field available on ThALES, we performed one more experiment at the FLEXX spectrometer using the 15 T magnet equipped with a dysprosium "booster" (setup 2), which has a much smaller sample volume, so the sample mass and, consequently, the count rate were reduced (which was partly compensated by choosing a larger value of k f = 1.5 Å −1 ). The resulting spectra, shown in Fig. 2 (b), had to be counted up to 2 h per data point. In this configuration, we could reach fields up to 16.5 T, following a shift in the energy of resonance C with approximately the same slope (g-factor) as that of the more intense resonance A. This leaves no doubt about the magnetic origin of the second peak.
In Fig. 3 (a), we summarize these new data (setups 1-3) together with the previously published results for the same field direction (setups 6, 10, 11) [43] in the form of a color map. The fitted positions of resonances A and C are shown with red diamonds and black pentagon symbols, respectively. One can see that the energy of both resonances shifts upwards with increasing field with approximately the same slope. The expected position of resonance B, according to ESR results [41], is shown with a dashed gray line. In the low-field region within Phase III, an arrow marks the previously reported feature at twice the AFM charge gap [62].

Dispersion of field-induced multipolar excitations
In the two bottom panels of Fig. 3, similar field dependencies for the same magnetic field direction B [110] are presented for two other high-symmetry points in the BZ. The data at R( 1 Fig. 3 (e) combine our previously published results for fields below 6 T (setups 6 and 10 from Refs. [62] and [43], respectively) with new measurements using setups 1 and 3 that extend to higher fields up to 14.5 T. The corresponding raw data for the three highest field values are also presented in Fig. 4 (a). One can see that the two excitations that were observed earlier in phase II [43] persist also in higher fields, following the same linear trends with an approximately equal slope. The data at M ( 1 2 1 2 0) in Fig. 3 (f) consist of the zero-field and 2.5 T spectra from Ref. [43] (setup 10) [59] and four new scans at higher fields between 8.5 and 14.5 T, which are shown separately in Fig. 4 (b). They indicate the presence of at least one excitation whose energy increases monotonically with the magnetic field.  . One can see that the more intense higherenergy excitation at the R point appears to be continuously connected to resonance A at the zone center, whereas the weaker low-energy resonance at R seems to cross this branch and re-emerge as resonance C at the Γ point. At the same time, both resonances approach each other along the Γ M line, so that only a single peak is observed at the M point. However, from our earlier results [43], we know that the low-energy peak at the R point emerges already at 1.7 T, simultaneously with phase II. On the other hand, resonance C at the Γ point is only visible above 14 T, which leaves an open question about the detailed evolution of this branch along the RΓ line with increasing field. Unfortunately, the two excitations cannot be clearly resolved close to their crossing point between R and Γ , hence it is unlikely that similar momentum-dependent maps at intermediate magnetic fields will help to answer this question conclusively.
Finally, it is useful to compare the dispersion of multipolar excitations in Fig. 5 (a) with similar momentum-resolved intensity maps measured in smaller magnetic fields of 2.5 and 5 T, which were presented earlier [43]. Across the whole field range, the bottom of the magnon band is located at the Γ point, and the maximum of the dispersion is reached at the M point. The total magnon bandwidth does not change and stays at approximately 1.2 meV in the whole field range, while the whole spectrum shifts upwards in energy with increasing field. At 2.5 T, there are two resonances at the R point, of which the lower one displays a steep upward dispersion along the RΓ line, ultimately merging with the more intense branch emanating from the zone center. As we can see now, these two branches actually cross each other, leading to the appearance of resonance C at the Γ point that can only be seen in high magnetic fields. This behavior is qualitatively consistent with the theoretical predictions that will be presented in Sec. III C.

C. Anisotropy with respect to the field direction
We proceed by demonstrating the dependence of the INS spectra on the direction of the applied magnetic field. Figure 6 shows the comparison of raw INS data measured in the vicinity of the Γ (110) point at the same or similar values of magnetic field, but for different field orientations. Additional spectra are also presented in Figs. S2 and S1 of the Supplemental Material [63]. At a relatively low field of 4 T [ Fig. 6 (a)], the spectra look nearly independent of the field direction, except for the difference in absolute intensity and the signal-to-background ratio, which can be explained by different experimental setups and by the different orientation of the crystal rod with respect to the beam [64]. At higher magnetic fields [ Fig. 6 (b-d)], the differences become more pronounced. This is seen not only in the different energy of the resonance as a function of field direction, but also in the appearance of additional spectral lines for B [001]. For this field orientation, one can see three peaks that shift as a function of field with a different slope, as shown in Fig. S2 in the Supplemental Material [63] and in Fig. 3 (d). The same three peaks are also reproduced near the equivalent zone-center point Γ (100), as shown in the bottom curve of Fig. 7 (b). This example clearly demonstrates that INS measurements in a rotating field offer additional information about collective modes in CeB 6 that have so far avoided direct observations in either INS or ESR experiments. Qualitative differences are also seen at other points in the BZ. For instance, the dispersion along the Γ M direction, measured using setup 4 [57], is depicted in Fig. 5 (b). An additional peak is observed at the M point, which appears to correspond to two excitations that split closer to the zone center.
The measurements presented so far were conducted using different experimental setups with different cryomagnets and at somewhat different temperatures, depending on the base temperature of the available sample environment. These temperatures were always much lower than the ordering temperature of phase II, and therefore one expects that these small temperature differences can be neglected. Nevertheless, for a reliable comparison of the data acquired with different setups, we had to ensure that the possible temperature dependence of the excitation energies is not detectable within our experimental accuracy. This is demonstrated in Fig. 7, where measurements at the base temperature of 0.5 and 1.7 K in setups 6 and 4, respectively, are compared with similar spectra measured at an elevated temperature of 3.3 or 3.2 K . For both field directions, B [110] and B [001], the spectrum shows no detectable changes with temperature. In Fig. 7 (b), the higher background of the higher-temperature curve is a consequence of non-magnetic scattering, as this spectrum was measured exactly at the allowed Bragg peak position, while the low-temperature data were measured with a small offset in Q. Nevertheless, the Fig. 6. Comparison of the raw INS spectra at the Γ (110) point (or very close to it, to avoid Bragg contamination), measured in magnetic fields applied along different crystallographic directions. The reciprocal-space vector, measurement temperature, magnetic field and its direction are indicated in the legends. Note that the data were measured using different spectrometers in different configurations, therefore the absolute intensity is not directly comparable. The curves are shifted vertically for clarity by an amount indicated on the vertical axis with horizontal bars of corresponding colors. Solid lines are empirical fits that serve as guides to the eyes. positions of all three magnetic peaks in these spectra are identical within the resolution of our experiment. This way, we can exclude the effect of small temperature variations on the outcome of our measurements.
At first glance, this conclusion may seem to be in contradiction with the recent ESR results by Semeno et al. [65][66][67]. In this experiment, the energy of the zone-center magnetic resonance has been measured for different directions of the magnetic field between B [110] and B [001], revealing a significant (∼10%) anisotropy in the g-factor with a strong temperature dependence, which was most pronounced for B [001]. The g-factor measured in ESR spectroscopy can be directly related to the slope of the linear field dependence of magnetic excitations at the Γ point, and we have shown [43] that the values obtained from ESR and INS measurements for the field direction B [110] are in perfect agreement with each other. On the other hand, for B [001], the ESR results suggest a 20% reduction of the g-factor between T = 1.8 and 3.2 K that can be clearly excluded based on our neutron-scattering data in Fig. 7 (b). This apparent contradiction is resolved by plotting the position of the ESR line for different temperatures and overlaying it on the phase diagram of CeB 6 for the same field direction [68], see Fig. 8. Due to the technical restrictions of the ESR technique, the measurements have to be performed in a resonator with a fixed frequency (in this case, 60 GHz or 0.25 meV), and the microwave absorption is measured as a function of magnetic field. The red rectangle in Fig. 8 marks the region in the B-T parameter space covered by the data in Refs. [65][66][67], while the blue squares mark the positions of the 60 GHz resonance at different temperatures. In this representation, it becomes obvious that the field scans in the low-temperature region cross the phase boundary between the AFM and AFQ phases, where the magnetic excitation spectrum becomes quasielastic [43] and gets broadened by the critical fluctuations of both order parameters. In the proximity to a phase transition, the dependence of the resonance energy on magnetic field loses its linearity, which prohibits a reliable determination of the g-factor from a single measurement at one fixed frequency. We can therefore conclude that the temperature dependence of the ESR spectrum reported in Ref. [65] results primarily from the proximity to phase III, rather than from a change in the g-factor within phase II, which we have shown to be temperature independent. Moreover, given that there are at least three magnetic resonances revealed in our INS measurements [ Fig. 3 (d)], which cannot be clearly resolved at small fields, the ESR g-factor should correspond to an average slope of all three resonances. Taking the g-factor value of ∼ 1.4 at the maximal temperature of 3.2 K, which should be more reliable because the AFM phase is fully suppressed at this temperature, we find it to be in reasonable quantitative agreement with the average slope of the three INS resonances. We hope that these considerations will motivate additional ESR experiments in higher magnetic fields, which can be then directly compared with the presented INS data for B [001].

A. Localized model for Ce 4 f 1 and multipolar HO
We will now turn to the theoretical interpretation of the presented results and the discussion of the corresponding theoretical model that can explain the observed dependence of multipolar excitations on the magnitude and direction of the magnetic field. The localized 4 f model for the HO in CeB 6 has been proposed [69] and investigated with respect to possible multipolar order [8]. Its basic features will be briefly summarized here as a basis for interpreting experimental findings. It starts from the observation that the CEF splitting of 4 f 1 Ce 3+ Γ 8 ground state quartet and excited Γ 7 state is very large and therefore the latter may be neglected at low temperatures [70][71][72]. In terms of free 4 f 1 (J = 5 2 ) states |M 〉 (|M | ≤ J) the quartet states |τσ〉 (τ = ±, σ =↑↓) are given by The fourfold degenerate ground-state manifold may be thought of being composed of two doublets of different Γ 6 and Γ 7 "orbital" symmetry, each being twofold Kramers degenerate. Therefore they can be presented by two types of pseudo-spins, namely τ for the two different orbitals and σ for the Kramers degeneracy. The complete set of operators with n = 1-15 and α, β = x, y, z is constructed from the Pauli matrices of pseudo-spins. They constitute a basis set for the fifteen multipolar moments, acting in the space of Γ 8 quartet states, which are irreducible combinations of the X n in O h cubic symmetry. In terms of these variables the model Hamiltonian describing the effective intersite coupling of Γ 8 multipoles may be written as The first term is a SU(4) "supersymmetric" nearest-neighbor intersite interaction on the simple-cubic Ce-sublattice (sites i, j) of strength D, which has no bias for any of the fifteen multipoles as primary order parameter. The second term describes the symmetry-breaking terms that favor quadrupolar or octupolar order depending on the size of the parameters = (ε Q , ε O ), which correspond to Γ + 5 -type quadrupoles and Γ − 2 -type octupoles, respectively. The ± denotes even/odd behavior under time reversal. The last term is the Zeeman energy in pseudospin representation, where H = B/µ 0 is the external magnetic field. The equivalence of multipole representation in terms of conventional total angular momentum operators J and the pseudospin expressions is presented in Table II.
The previous conclusion from critical field anisotropy [8], field-induced neutron diffraction [21], NMR results [73], and resonant x-ray scattering [23] indicate that HO may be well described as a primary AFQ Γ + 5 order with a propagation vector Q = (π, π, π) and a secondary strongly field-induced octupolar Γ − 2 order parameter at the same wave vector. This may be achieved by choosing = (0.5, 0.5) [28] and should be considered as an appropriate set for CeB 6 [71].

B. RPA calculation of magnetic excitations in the multipolar-ordered phase
For the model in Eq. (3), the magnetic excitation spectrum of the multipolar-ordered phase as measured by INS may be calculated with complementary generalized Holstein-Primakoff approach [28] or multipolar response function formalism in the RPA approach [29,30]; the latter will be used here. As the first step, one has to calculate the effective molecular fields 〈X n A 〉 ± 〈X n B 〉 (uniform and staggered) of each multipole basis operator where s = A, B denote the simple-cubic sublattices of the antiferro-type HO defined by ordering vector Q. The CEF states are mixed by the molecular fields into new eigenstates |νs〉 i with energies E s ν at every sublattice site (s, i). In terms of standard basis operators a s νµ (i)= |νs〉 i 〈µs| i (ν, µ = 1, 4) for the molecular field states degeneracy O h rep., multipole Stevens notation pseudo-spin form Table II. Table of 7 of the 15 multipoles of Γ 8 quartet in different representations: Stevens notation using total angular momentum components J α , α = x, y, z, or pseudospin components σ α , τ α . In the first three rows we used the multipole vector η = (−τ z σ x + 3τ x σ x , − τ z σ y − 3τ x σ y , 2τ z σ z ). The components of the total angular momentum are the linear combination J α = n λ n α X n where coefficients λ n α can be read off in the last column. In the last row, symmetrization (summation over all permutations of x, y, z) is denoted by a bar. the mean field approximation to Eq. (3) is then given by where M is a 15-component vector of matrix elements for the multipole operators defined by M ns νµ = 〈ν, s|X n is |µ, s〉 and the inter-sublattice multipole 15 × 15 diagonal interaction matrix is D AB (q) = D BA (q) = −2z Dγ q Λ. Here γ q = 1 3 (cos q x + cos q y + cos q z ), z = 6, and Λ(n, n ) = Λ(n)δ n,n with Λ(5) The collective response of all 15 Γ 8 multipoles for the two sublattices is then described by the 30×30 RPA susceptibility matrix From its elements we construct the physical-moment (J) cartesian 3 × 3 susceptibility matrix according to where λ n α are the coefficients of J α in the pseudo-spin representation (Table II). Then the dynamical dipolar structure function seen in INS may be written as , ω; B), (8) which depends parametrically on field strength and direction.
Here Q = q + K is the total momentum transfer in INS with K denoting a reciprocal lattice vector andQ = Q/|Q|. The prefactor under the sum projects to the contributions with momentum transfer perpendicular to the dipolar moment J.
In the following the structure function will be computed for various model parameters, field strengths and continuous rotations of field direction and compared to the experimental INS scattering results.

C. Dispersion of multipolar modes in phase II
First we discuss the dispersion of multipolar modes for a fixed magnetic field B. In the model calculations, we only include the AFQ order of phase II, the effect of the complicated dipolar magnetic order (phase III) [21,71] at a different wave vector is not treated here. An attempt to this purpose has been made in [74]. The model parameters are given by = (0.5, 0.5) [28,30], corresponding to the composite degenerate Γ + 5 quadrupole (O yz , O z x , O x y ) and Γ − 2 (T x yz ) octupole order. It should be noted that the precise values of ε Q , ε O are not known. The accidental degeneracy ε Q = ε O is assumed to minimize the number of parameters. This model is known to reproduce the magnetic-fieldtemperature phase diagram of CeB 6 and its field anisotropies well, including the fluctuation effects beyond the mean-field approximation [75] as well as the essential NMR results [73]. In particular, ε Q = ε O can explain the strongly convex shape of induced octupole T x yz (B) field dependence observed directly in RXS experiments [76].
Furthermore k B T 0 = 2zD = 0.41 meV or T 0 = 4.74 K defines the interaction energy and temperature scale. Following the notation of Ref. [30], we can also introduce the dimensionless field strength h = h/T 0 with h = (7/6)g J µ B H/T 0 (assuming k B ≡ 1), which can be expressed as h = 0.672B[T]/T 0 [K] using physical units for B and T 0 . Note that the mean-field transition temperature of the model is T mf Q (h = 0) = (1 + ε Q )T 0 = 7.12 K, which is considerably larger than the experimental transition temperature T mf exp = 3.2 K, because the latter is strongly suppressed by fluctuations [75].
The calculated dynamical RPA dipolar structure function S(Q, ω) for CeB 6 is shown in Fig. 9. The polygonal path in momentum space and the strength and magnitude of the magnetic field used in the calculation were chosen to enable comparison with the experimental data. The dispersion of multipolar excitations is in excellent agreement with the results of the alternative HP approach, including the intensities of modes as a function of Q [28,30]. Because there are two sublattices and locally three excited CEF states, in the AFQ/AFO-type ordered phase six excitation branches will ap- pear prominently at low temperature that are mostly visible, e.g. in Fig. 9 (b-d) (for h = 0 there are additional degeneracies). Roughly the six branches can be arranged in two groups: The low-energy branches with a strong field dependence are stabilized directly by the Zeeman term, whereas the remaining high-energy branches are mainly stabilized by the strongly field-induced octupolar molecular field [28].
When temperature approaches T mf Q (h ) from below at a constant field, the high-energy modes are expected to collapse due to the reduced octupolar order, while the lowenergy modes should be less affected. Complementary, when temperature is kept constant much below T mf Q and the field is reduced to zero, the high-energy modes change little, while the low-energy modes are shifted downwards. Due to the threefold degenerate Γ + 5 order, a Goldstone mode then appears for h = 0 at the Γ point, see Fig. 9 (a). At finite field, the degeneracy is lifted by choosing a linear combination of the three Γ + 5 quadrupoles as the ground state (see Sec. III D 2), leading to a Γ -point anisotropy gap in the excitation spectrum in Figs. 9 (b-d). It should be mentioned that the relative field independence of the higher-energy modes is a consequence of the accidental degeneracy ε Q = ε O assumed in the model [30] which is consistent with the RXS results [76].
The color plots in Fig. 9 show not only the dispersion, but also the expected intensity of the multipolar excitations that can be directly compared with the INS experiments presented in Sec. II B 2. The INS data in Fig. 5 (a), measured along the R -Γ -M path for B [110], should be directly compared with the left half of Fig. 9 (d), whereas the experimental color map in Fig. 5 (a), measured along the Γ -M direction for B [001], should be compared with the corresponding section of Fig. 9 (b). As mentioned in the introduction, we expect no quantitative agreement with experiment away from the Γ point, as our theoretical model neglects important interactions beyond the nearest neighbors, yet even such an oversimplified approach reveals some qualitative similarities with the measured data. For instance, the crossing of the two low-energy branches along the Γ -M direction for B [001] and along the R -Γ direction for B [110] are quite well captured by the model.
Note that under an external field, the intensity distribution no longer follows the cubic symmetry of the lattice. The dynamic structure factor depends on the relative orientation of the Q and B vectors, as illustrated in Figs. S4 and S3 in the Supplemental Material [63]. For instance, the intensity of one of the high-energy modes at Q = Γ (100) fully vanishes if the magnetic field is applied parallel to the Q vector, B [100], but is finite for equal but orthogonal field B [001] ⊥ Q (Fig. S4). Similar changes in intensity are seen near Q = Γ (110) for magnetic fields along [110] ⊥ Q or [110] Q (Fig. S3).

D. Dependence of multipolar excitations on the field strength and direction
We will now consider the magnetic-field dependence of mode energies for a fixed momentum transfer Q. To minimize the number of free parameters, we will restrict our consideration to the high-symmetry points of the BZ that were investigated experimentally. As pointed out in the introduction, the benefit of this approach is that for a fixed Q corresponding to a symmetry point, the various interactions enter only in specific sums and not their individual values which should be much less sensitive to the unavoidable inaccuracies in the description of the multipolar interaction Hamiltonian. For a fixed Q, the momentum-dependent part of the model can result in an energy offset of the field dependence, which we expect to vanish at the zone center (Γ point) at least for the low-energy excitation branches. Looking at the multipolar excitations in field space (magnitude and direction of B) for a fixed momentum is a change of viewpoint that represents a complementary approach as compared to the previous theoretical [28,30] and experimental [34,43] investigations.

Variation of mode energies with field strength
First, we consider the continuous dependence of mode energies at the high-symmetry points of the BZ on field strength, which is illustrated in Fig. 10 . This corresponds to a rotation of the field direction in the plane orthogonal to the Q vector. In addition, a similar field dependence for B Q (which was not investigated experimentally) is presented in Fig. S6 (d) of the Supplemental Material [63].
There are very pronounced differences in the mode energies for different field directions, related to the behavior of the O x y quadrupoles and the Γ − 2 octupoles, that form primary and secondary (field-induced) order parameters of phase II, with an applied field. The O x y and equivalent O yz and O z x quadrupoles are selected from the Γ + 5 manifold as the primary AFQ order parameter. Their threefold degeneracy (Table II) enables a free rotation of the charge density of the quadrupole such that its lobes always point in the directions orthogonal to the field, as shown in Fig. 11. For a magnetic field pointing along the direction defined by the unit vector H/|H| = (αβγ), the primary order parameter can be written as the linear combination αO yz +βO z x +γO x y [71]. In contrast, the field-induced octupole T x yz is nondegenerate and therefore fixed with respect to the crystal axes when H rotates.
As mentioned before, the high-energy modes ω 4T 0 ≈ 1.6 meV are stabilized by the octupolar molecular field and are hardly influenced by the applied magnetic field. They can be seen as horizontal, approximately field-independent lines in Fig. 10. On the other hand, the low-energy modes are stabilized directly by the Zeeman term and show an approximately linear increase in frequency and splitting with the field, with the slope given by the effective g-factor that can be anisotropic with respect to the field direction [65][66][67]. At sufficiently high fields, these modes intersect with the high-energy modes, which leads to a complex hybridization pattern with multiple avoided mode crossings. This can lead to drastic changes in the slope and a redistribution of intensity between different modes above a certain field threshold, possibly leading to the appearance of an additional ESR line [41] and to the second higher-energy peak in the INS spectrum described in Sec. II B.
To simulate how the specific choice of the Q vector at which the Γ -point spectra are measured in an INS experiment influences the relative intensity of the multipolar modes, in Figs. S5 and S6 in the Supplemental Material [63] we compare the calculated field dependencies for B 2 0), that were investigated experimentally. In spite of the limited INS data for these points, it is clear that the linear field dependence at the R point in the experimental data ( Fig. 3) extends to much higher fields and does not exhibit the pronounced flattening of the strongest excitation near 8 T, suggested in the calculations. We attribute this discrepancy to the limitations of our model, which is not expected to be quantitatively accurate at wave vectors other than the zone center.

Multipolar excitations under rotation of the field direction
Another result of our new approach to the multipolar excitations in CeB 6 concerns the dependence on continuous field rotation and is presented in the polar graphs of Fig. 12 for h = 2 (B ≈ 14 T). In these figures, the radial coordi- Fig. 11. Sketch of the quadrupole (left) and octupole (right) order parameter in a rotating field that points along the unit vector (αβγ) obtained from the real-space equivalent tesseral harmonics of the Stevens operators in Table II. The distance from center to the lobes is proportional to the charge/spin-density of the even Γ + 5 quadrupole or odd Γ − 2 octupole. The total charge or spin moment is zero. Therefore equal regions with ± charge or (N/S) spin densities (different colors for each) appear. When the field unit vector (αβγ) rotates, the quadrupole charge density rotates with it such that it remains in the plane perpendicular to H. This is due to the 3-fold degeneracy of Γ + 5 . In contrast, the non-degenerate Γ − 2 octupole spin density remains fixed with respect to the crystal axes when H rotates. As a function of field strength H, the primary quadrupole order parameter is little affected, while the induced octupole shrinks to zero for H → 0. For each field direction given by the unit vector (αβγ), the AFQ order parameter αO yz + βO z x + γO x y is selected from the Γ + 5 manifold [71] (see Fig. 11), which is automatically included in the mean-field ground state of the model given by Eq. (3). The resulting field-angular mode anisotropies show various signatures that are worthwhile to mention and to look for in future experiments. First of all, the anisotropy of the low-energy modes that exhibit approximately linear field dependence at small fields represents the effective g-factor anisotropy that can be directly measured both in ESR (for the Γ point) [65][66][67] and in INS (for any Q vector). The highenergy mode frequency (large polar radii) is less sensitive to the magnetic field and therefore changes (expands) slower with increasing field strength than the low-energy modes (small polar radii). As a result, an increase in field strength leads to a hybridization of the low-and high-energy modes, which generally increases the complexity and anisotropic appearance of the field-angular variation of mode frequencies. We also note that the anisotropy patterns at the Γ and R points are quite similar. This is expected since both points have full cubic symmetry. Interestingly, the relative intensity of the low-energy and high-energy modes is interchanged when going from Γ to R and vice versa, in agreement with the experimental result in Fig. 5 (a). At the X and M points, the low-energy mode frequencies are larger (inner radii increased), but again we observe some intensity exchange, now among the low-energy modes themselves.
Altogether these polar field-angle anisotropy plots of multipolar excitations in CeB 6 at high-symmetry points represent a large set of information on mode positions and intensities that may be very useful for comparing with future experimental results and should give guidance where to look for the modes with largest intensity. If the model so far accepted for CeB 6 with = (0.5, 0.5) is reasonable, some features of the field-anisotropy plots described above should be identified in INS experiment with quasi-continuous field rotation. Such a comparison becomes possible using the data presented in Sec. II C.

E. Comparison to the experimental results
A reliable comparison with experiment is so far only possible at the Γ point, where detailed field-angular dependences of the INS spectra were obtained (see, for instance, Figs. 1, 2, 3, 6, 7). In Fig. 13, the corresponding experimental peak positions for the fields of 4, 7, 10, and 14.5 T are plotted on top of the calculated field-angular polar maps for the same magnetic fields. Every segment of the plot covers the same irreducible 90 • range of field directions between [001] and [110] in the plane orthogonal to [110]. The experimental data points were obtained either by directly fitting the peak positions from measurements at the corresponding field values (filled symbols) or via linear interpolation of the field dependences (open symbols). Different symbol shapes for different modes are consistent with those in Fig. 3.
We see a remarkably good agreement of the calculation with the field-directional anisotropy of the most intense low-energy mode (resonance A) with a quasi-linear field dependence, which has been followed experimentally over a large field range for all four high-symmetry directions of the magnetic field: B [001], [112], [111], and [110]. At high magnetic fields of the order of 10 T, we observe a very considerable anisotropy of ∼ 60% in the effective g-factor between its minimal and maximal values reached for the [001] and [111] field directions, respectively. For lower fields, however, the anisotropy is reduced, which is a direct consequence of the initially nonlinear field dependence of the low-energy mode that can be seen in Fig. 10 (d). Indeed, the g-factor anisotropy reported from previous ESR measurements [65][66][67], which were performed at magnetic fields of ∼ 3 T, is several times smaller and constitutes a relative change of less than 10% between the [100] and [111] field directions. However, the overall shape of the anisotropic g-factor dependence turns out to be the same in ESR and INS measurements, if one considers only those ESR data that were measured at an elevated temperature away from the transition line to the AFM phase (see Fig. 8 and the discussion at the end of Sec. II C). The same anisotropic field-angular dependence could be well described by a model proposed in Refs. [77][78][79] by Schlottmann.
Another experimental confirmation of the possibly nonlinear field dependence of the low-energy mode in Fig. 10 (d) for the field direction B [001] follows from the experimental field dependence plotted in Fig. 3 (d). This is the only field direction, for which an extrapolation of the high-field part of the field dependence (measured in phase II, far from the AFM phase) for the position of the most intense lowenergy peak results in a small but significant offset on the energy axis, in agreement with the theoretical model. For other field directions, such as B [110] in Fig. 10 (a), the deviations from a linear dependence first start at much higher fields due to the hybridization with high-energy modes. Figure 13 also shows the positions of the new resonance peaks that have been so far measured only for B [110] and B [001] in high magnetic fields. Their energies fall in the range where the theoretical model predicts multiple modes resulting from the hybridization of high-and low-energy excitations. With the available data, it is not possible to assign one of these modes uniquely to the experimentally observed resonances. However, the field dependence in Fig. 10 (a) suggests that the high-energy mode (resonance C) may correspond to the upper branch that gains intensity beyond the avoided mode crossing for B 10 T and hω 2 meV. This field and energy range approximately corresponds to the region where resonance C has been seen in our experiment. The predicted energy splitting between the two modes of ∼ 1 meV also roughly matches with the energy difference between the INS resonances A and C.
On the other hand, the model predicts no additional resonances at the Γ point below the energy of resonance A, which means that it is unable to explain the appearance of the resonance B in ESR data at fields above 12 T [41]. Theoretical considerations about the possible origin of this high-field resonance were proposed earlier by Schlottmann [77][78][79], but the reason why it cannot be seen in INS measurements remains unclear. The dispersion of resonance B, as well as its relationship to the multipolar excitations described in our work, also remain to be clarified.

IV. SUMMARY AND DISCUSSION
In summary, we have proposed an alternative approach to study magnetic excitations in a correlated-electron system characterized by an interplay of charge and orbital degrees of freedom, based on the measurements of INS spectra at fixed points in momentum space as a function of magnetic field and its direction. This approach is complementary to the conventional analysis of INS data, which primarily considers the dispersion of magnetic excitations as a function of momentum for a single magnetic-field direction. Measurements of the field-directional anisotropy are a well-established technique in many physical probes, such as thermodynamic or transport measurements, and it is therefore natural to extend this approach to the local spectroscopic probes, such as neutron spectroscopy. In complex systems characterized by strong spin-orbit coupling and exhibiting multipolar order parameters, INS measurements of the field-angular anisotropy of magnetic excitations can help in disentangling the "local" physics related to the fielddependent molecular field and the Zeeman splitting of the single-ion eigenstates from the effects of magnetic interactions that lead to the formation of dispersive collective modes with a pronounced momentum dependence.
We have applied the proposed method to the analysis of INS data on the long-studied material CeB 6 , which exhibits complex low-temperature ordering phenomena and rich magnetic excitation spectra that have so far avoided a complete theoretical description in spite of many years of experimental and theoretical effort. A comparison of our experimental results presented in Sec. II with the model calculations described in Sec. III reveals reasonable agreement of the excitation spectra in the vicinity of the zone center (Γ point). In particular, the available model can reproduce the experimental g-factor anisotropy and can qualitatively explain the appearance of the new high-energy mode in the excitation spectrum that can only be observed at high magnetic fields. On the other hand, the model is much less accurate in describing the dispersion of multipolar excitations away from the zone center, as it neglects significant long-range RKKY couplings between the next-nearest and next-next-nearest neighbors. Potentially, combining INS data in momentum and field space may provide sufficient information to fit multiple interaction terms in a more realistic Hamiltonian of CeB 6 in a similar way, as it is routinely done for conventional dipolar magnets.
The vector field effectively adds three new dimensions to the parameter space of the INS measurement, making it exceedingly difficult to perform the full mapping of the parameter space, as it would require an unreasonably large amount of neutron beam time. This is perhaps one of the reasons why such an approach has not been applied to any other material until now, to the best of our knowledge. Another reason is more technical, related to the development of stable and reliable mechanical rotators that would be compatible with high-field magnets used in INS experiments [80]. In future, with the advent of new high-flux neutron sources and the development of suitable sample environment for the stable mechanical rotation of single crystals around a specific scattering vector, measurements similar to those presented here for CeB 6 can be accomplished much more efficiently and for a broad range of materials. The subsequent comparison of such measurements with the corresponding theoretical model of multipolar excitations is able to provide the most direct information about inter-multipolar interactions that stabilize magnetically hidden order and to probe the HO symmetry across the magnetic phase diagram. We hope that our present work will guide future developments in this direction.    Fig. S5. Comparison of the field dependencies for the dynamic structure function at the Γ point in a magnetic field applied along the [001] cubic axis, calculated for different reciprocal-space points as indicated in the corresponding panels. In panels (a, b), the field is assumed to be perpendicular to the momentum-transfer vector, B ⊥ Q, in panel (d) the two vectors are parallel, B Q, whereas in panel (c) the angle between the Q vector and the field direction is ∼ 46 • . Fig. S6. Comparison of the field dependencies for the dynamic structure function at the Γ point in a magnetic field applied along the [110] axis, calculated for different reciprocal-space points as indicated in the corresponding panels. In panels (a-c), the field is assumed to be perpendicular to the momentum-transfer vector, B ⊥ Q, whereas in panel (d) the two vectors are parallel, B Q.