Magnon interactions in the frustrated pyrochlore ferromagnet Yb$_2$Ti$_2$O$_7$

TL;DR

This paper investigates magnon interactions in the frustrated pyrochlore Yb$_2$Ti$_2$O$_7$, revealing their significant role in extending spin-wave theory and explaining experimental observations of its complex excitation spectrum.

Contribution

It demonstrates how magnon interactions influence phase boundaries and excitation spectra, improving the understanding of Yb$_2$Ti$_2$O$_7$'s magnetic behavior across different fields.

Findings

Magnon interactions shift the phase boundary near the ferromagnet and $ ext{Gamma}_5$ phase.

Corrections to critical fields align theory with experiments.

Partial reproduction of low-field broadening in neutron scattering data.

Abstract

The frustrated rare-earth pyrochlore Yb$_2$Ti$_2$O$_7$ is remarkable among magnetic materials: despite a ferromagnetically ordered ground state it exhibits a broad, nearly gapless, continuum of excitations. This broad continuum connects smoothly to the sharp one-magnon excitations expected, and indeed observed, at high magnetic fields, raising the question: how does this picture of sharp magnons break down as the field is lowered? In this paper, we consider the effects of magnon interactions in Yb$_2$Ti$_2$O$_7$, showing that their inclusion greatly extends the reach of spin-wave theory. First, we show that magnon interactions shift the phase boundary between the (splayed) ferromagnet (SFM) and the antiferromagnetic $\Gamma_5$ phase so that Yb$_2$Ti$_2$O$_7$ lies very close to it. Next, we show how the high-field limit connects to lower fields; this includes corrections to the critical fields for the $[111]$ and $[1\bar{1}0]$ directions, bringing them closer to the observed experimental values, as well as accounting for the departures from linear spin wave theory that appear in $[001]$ applied fields below 3 T [Thompson et al., Phys. Rev. Lett. 119, 057203 (2017)]. Turning to low-fields, though the extent of the experimentally observed broadening is not quite reproduced, we find a rough correspondence between non-linear spin-wave theory and inelastic neutron scattering data on both a single-crystal sample, as well as on a powder sample [Pe\c{c}anha-Antonio et al., Phys. Rev. B, 96, 214415 (2017)]. We conclude with an outlook on implications for future experimental and theoretical work on Yb$_2$Ti$_2$O$_7$ and related materials, highlighting the importance of proximity to the splayed ferromagnet-$\Gamma_5$ phase boundary and its potential role in intrinsic or extrinsic explanations of the low-field physics of Yb$_2$Ti$_2$O$_7$.

Equations

$$H$$

$$J_{\pm \pm }\left({\gamma }_{ij}{S}_r^{+}{S}_{r^{\prime }}^{+}+h.c.\right)+J_{z\pm }\left({\zeta }_{ij}\left[S_r^z{S}_{r^{\prime }}^{+}+S_r^{+}{S}_{r^{\prime }}^z\right]+h.c.\right)]$$

$$-\sum _r\mathbf{B}\cdot {\mathbf{\mu }}_r,$$

$${\mathbf{\mu }}_r\equiv -{\mu }_B\left[g_{\pm }\left(S_r^x{\hat{\mathbf{x}}}_r+S_r^y{\hat{\mathbf{y}}}_r\right)+g_z{S}_r^z{\hat{\mathbf{z}}}_r\right].$$

$${\gamma }_x=-{\zeta }_x^{\ast }$$

$$J_{zz}$$

$$J_{\pm \pm }$$

$$g_z$$

$$H\equiv \frac{1}{2}\sum _{\mathbf{r}{\mathbf{r}}^{\prime }}\sum _{\alpha {\alpha }^{\prime }}{\mathbf{S}}_{\mathbf{r}\alpha }^{⊺}{\mathbf{J}}_{\mathbf{r}-{\mathbf{r}}^{\prime },\alpha {\alpha }^{\prime }}{\mathbf{S}}_{{\mathbf{r}}^{\prime }{\alpha }^{\prime }},$$

$${\mathbf{S}}_{\mathbf{r}\alpha }\equiv$$

$$+\left(S-n_{\mathbf{r}\alpha }\right){\hat{\mathbf{e}}}_{\alpha ,0},$$

$${\mathcal{J}}_{\mathbf{\delta },\alpha {\alpha }^{\prime }}^{\mu {\mu }^{\prime }}\equiv {\hat{\mathbf{e}}}_{\alpha ,\mu }^{⊺}{\mathbf{J}}_{\mathbf{\delta },\alpha {\alpha }^{\prime }}{\hat{\mathbf{e}}}_{{\alpha }^{\prime },{\mu }^{\prime }},$$

$$H\approx NS(S+1){ϵ}_{\mathrm{cl}}+S{H}_2+S^{1/2}{H}_3+S^0{H}_4,$$

$${ϵ}_{\mathrm{cl}}\equiv \frac{1}{2N_s}\sum _{\alpha {\alpha }^{\prime }}\sum _{\mathbf{\delta }}{\mathcal{J}}_{\mathbf{\delta },\alpha {\alpha }^{\prime }}^{00},$$

$$H_2$$

$$H_3$$

$$H_4$$

$$A_{\mathbf{k}}^{\alpha \beta }={\mathcal{J}}_{\mathbf{k},\alpha \beta }^{+-}-{\delta }_{\alpha \beta }\sum _{\mu }{\mathcal{J}}_{0,\alpha \mu }^{00},$$

$$B_{\mathbf{k}}^{\alpha \beta }={\mathcal{J}}_{\mathbf{k},\alpha \beta }^{++},$$

$$T_{\mathbf{k}{\mathbf{k}}^{\prime }}^{\alpha \beta \mu }=-\left[{\delta }_{\alpha \mu }{\mathcal{J}}_{{\mathbf{k}}^{\prime },\beta \alpha }^{+0}+{\delta }_{\beta \mu }{\mathcal{J}}_{\mathbf{k},\alpha \beta }^{+0}\right],$$

$$V_{\mathbf{k}{\mathbf{k}}^{\prime }[\mathbf{q}]}^{\alpha \beta \mu \nu }=\left({\delta }_{\alpha \mu }{\delta }_{\beta \nu }{\mathcal{J}}_{\mathbf{k}-{\mathbf{k}}^{\prime }+\mathbf{q},\alpha \beta }^{00}+{\delta }_{\alpha \nu }{\delta }_{\beta \mu }{\mathcal{J}}_{\mathbf{q},\alpha \beta }^{00}\right)-\left({\delta }_{\mu \nu }{\delta }_{\mu \beta }{\mathcal{J}}_{\mathbf{k}+\mathbf{q},\alpha \nu }^{+-}+{\delta }_{\alpha \beta }{\delta }_{\alpha \mu }{\mathcal{J}}_{\mathbf{k},\alpha \nu }^{+-}\right),$$

$$D_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}^{\alpha \beta \mu \nu }=-\frac{3}{4}\left({\delta }_{\alpha \mu }{\delta }_{\alpha \nu }{\mathcal{J}}_{{\mathbf{k}}^{\prime },\beta \alpha }^{++}+{\delta }_{\mu \beta }{\delta }_{\nu \beta }{\mathcal{J}}_{\mathbf{k},\alpha \beta }^{++}\right),$$

$${\mathcal{S}}_{\nu {\nu }^{\prime }}(\mathbf{q},\omega )\equiv \int d\omega e^{i\omega t}⟨{\mu }_{-\mathbf{q}}^{\nu }(t){\mu }_{\mathbf{q}}^{{\nu }^{\prime }}⟩,$$

$$\mathcal{S}(\mathbf{q},\omega )\equiv \sum _{\nu {\nu }^{\prime }}\left({\delta }_{\nu {\nu }^{\prime }}-{\hat{q}}_{\nu }{\hat{q}}_{{\nu }^{\prime }}\right){\mathcal{S}}_{\nu {\nu }^{\prime }}(\mathbf{q},\omega ),$$

$$I(\mathbf{q},\omega )\propto F(\mathbf{q}{)}^2\mathcal{S}(\mathbf{q},\omega ),$$

$${\mathbf{G}}^R(\mathbf{k},\omega )$$

$$\bm{M}_{\bm{k}}=\left(\begin{array}[]{cc}\bm{A}_{\bm{k}}&\bm{B}_{\bm{k}}\\ {\bm{B}}^{*}_{-\bm{k}}&{\bm{A}}^{*}_{-\bm{k}}\end{array}\right),$$

$$J$$

$$K$$

$${\hat{\mathbf{z}}}_1$$

$${\hat{\mathbf{z}}}_2$$

$${\hat{\mathbf{z}}}_3$$

$${\hat{\mathbf{z}}}_4$$

$$\left(\begin{array}[]{ccc}J+K&+D/\sqrt{2}&+D/\sqrt{2}\\ -D/\sqrt{2}&J&\Gamma\\ -D/\sqrt{2}&\Gamma&J\end{array}\right),$$

$$\left(\begin{array}[]{c}J\\ K\\ \Gamma\\ D\end{array}\right)=\frac{1}{3}\left(\begin{array}[]{cccc}-1&+4&+2&-2\sqrt{2}\\ +2&-8&+2&-2\sqrt{2}\\ -1&-2&-4&-2\sqrt{2}\\ -\sqrt{2}&-2\sqrt{2}&+2\sqrt{2}&+2\end{array}\right)\left(\begin{array}[]{c}J_{zz}\\ J_{\pm}\\ J_{\pm\pm}\\ J_{z\pm}\end{array}\right),$$

Full text

I Introduction

Of the hundreds of known three dimensional and (anisotropic) two dimensional magnetic materials, the overwhelming majority order magnetically at low temperatures and have magnetic excitations that may be understood as the quasiparticle modes corresponding to harmonic fluctuations of the underlying magnetic structure. One of the grand endeavours of condensed matter physics has been to identify materials where quantum fluctuations cause departures from this canonical picture. A particularly fruitful strategy to find such systems has been to focus on materials where the magnetic interactions are frustrated Lacroix et al. (2011), such as through geometrical frustration on lattices of corner-sharing triangles or through frustration arising from intrinsically anisotropic exchange couplings.

A noteworthy example is found in Yb2Ti2O7: a rare earth magnet with effective spin one-half ytterbium ions on the pyrochlore lattice of corner-sharing tetrahedra Gardner et al. (2010); Rau and Gingras (2019); Hallas et al. (2018); Blöte et al. (1969); Hodges et al. (2001a, b, 2002); Yasui et al. (2003); Bonville et al. (2004); Gardner et al. (2004); Malkin et al. (2004); Cao et al. (2009a, b); Ross et al. (2009); Thompson et al. (2011a, b); Onoda and Tanaka (2011); Onoda (2011); Yaouanc et al. (2011); Ross et al. (2011a); Chang et al. (2012); Hayre et al. (2013); Applegate et al. (2012); Ross et al. (2011b, 2012); D’Ortenzio et al. (2013); Chang et al. (2014); Lhotel et al. (2014); Pan et al. (2014); Jaubert et al. (2015); Robert et al. (2015); Gaudet et al. (2015); Pan et al. (2016); Bhattacharjee et al. (2016); Tokiwa et al. (2016); Gaudet et al. (2016); Yaouanc et al. (2016); Yan et al. (2017); Peçanha-Antonio et al. (2017); Kermarrec et al. (2017); Thompson et al. (2017); Arpino et al. (2017); Mostaed et al. (2017); Li et al. (2018); Shafieizadeh et al. (2018). The magnetic ground state of Yb2Ti2O7 has been directly found to be ferromagnetic with a net spontaneous moment along a cubic lattice direction Yasui et al. (2003); Chang et al. (2012); Yaouanc et al. (2016); Kermarrec et al. (2017), further supported by the existence of symmetry-breaking phase transitions in $[1\bar{1}0]$ Ross et al. (2011b) and $[111]$ magnetic fields Scheie et al. (2017) and by the absence of such a transition in a $[001]$ magnetic field Thompson et al. (2017). The finite temperature transition into the low temperature magnetically ordered state is known to be sensitive to the details of sample synthesis Blöte et al. (1969); Hodges et al. (2001a); Bonville et al. (2004); Yasui et al. (2003); Yaouanc et al. (2011); Ross et al. (2012); Lhotel et al. (2014); Gaudet et al. (2015, 2016); Mostaed et al. (2017) but samples have been made in which a sharp transition is observed in the range $200-260$ mK Chang et al. (2012); Thompson et al. (2017); Scheie et al. (2017).

In striking contrast to conventional magnets, the magnetic excitations in Yb2Ti2O7 in zero applied magnetic field are dominated by an apparent broad continuum of intensity. This behavior is reminiscent of the low-energy excitations observed in the (disordered) kagome antiferromagnet Herbertsmithite Han et al. (2012) and the higher-energy features seen in the (ordered) honeycomb antiferromagnet $\alpha$-RuCl3, which both exhibit anomalously broad excitation spectra Banerjee et al. (2016, 2017, 2018). However, unlike these materials, in Yb2Ti2O7 the exchange scale is small and the magnetic moments relatively large, so the magnetic ground state is tunable using easily accessible laboratory magnetic fields. This allows access to a high field regime where sharp one-magnon excitations are observed Ross et al. (2011b); Thompson et al. (2017).

In a single crystal with a sharp transition at $\sim 214\ {\rm mK}$, the magnetic excitations have been studied in magnetic fields along $[001]$ from high-fields ($9\ {\rm T}$) down to zero field Thompson et al. (2017). This experiment revealed the presence of sharp one-magnon excitations, well separated from the two-magnon excitations, that are well-described by linear spin-wave theory (LSWT) down to magnetic fields of about $3\ {\rm T}$. As the magnetic field was lowered further, the one-magnon dispersions were significantly renormalized from the LSWT predictions and level repulsion from the two magnon continuum was inferred. At low fields, as the two magnon continuum begins to overlap with the one magnon states one finds pronounced broadening effects. At zero field, the scattering intensity over large areas of momentum space was found to be a broad continuum within the resolution of the instrument.

These experimental results present two clear puzzles. Foremost is the question of how the magnetically ordered ground state and the continuum of excitations can be understood in a unified way that also smoothly recovers the conventional behavior observed at high-fields. Secondly, one would like to understand the origin of the sample dependence of the magnetic transition and whether this is related to the exotic zero field behavior. In this paper, we tackle the first of these questions by looking at the leading corrections to spin wave theory beyond the non-interacting theory.

The observed departures of the inelastic intensity from LSWT at intermediate fields are of a nature that one might expect to be described by extending LSWT to include the effects of magnon interactions. It is well understood that, generically, these interactions include terms that violate magnon number conservation - including three-magnon couplings Zhitomirsky and Chernyshev (2013); Chubukov et al. (1994); Chernyshev and Zhitomirsky (2009); Mourigal et al. (2013). These couplings induce spontaneous decay processes that mix the one-magnon and the two-magnon sectors. Such decays are kinematically forbidden when there is a separation in energy between the one- and two-magnon states (as is the case at high fields), only becoming operative at low fields when there is overlap between the two in energy-momentum space. Their effects have been extensively studied theoretically in various models of low-dimensional frustrated magnets where quantum fluctuations are a priori expected to be significant Zhitomirsky and Chernyshev (2013); Chubukov et al. (1994); Chernyshev and Zhitomirsky (2009); Mourigal et al. (2013); Plumb et al. (2016); Hwang and Kim (2016); Verresen et al. (2018, ). Experimentally, such effects have also been observed in several different magnetic materials using inelastic neutron scattering Masuda et al. (2010); Oh et al. (2013); Robinson et al. (2014); Dalla Piazza et al. (2015); Sandilands et al. (2015); Hong et al. (2017). In Yb2Ti2O7, the experimental evidence is clear that the onset of the broadening of one-magnon intensity coincides with the merging of the two-magnon continuum with the higher energy one-magnon branches Thompson et al. (2017). In addition, the observed broadening proceeds to lower energies as the field is further lowered and the multi-magnon states begin to significantly overlap the remaining sharp one-magnon modes.

As much as one might hope that non-linear spin-wave theory (NLSWT) can capture the intermediate field inelastic properties of Yb2Ti2O7, one also naively expects it to fail to describe the zero field features. This is because the computed magnon spectrum about the experimentally determined ferromagnetic state is gapped; the lowest energy multi-magnon states then lie at a finite threshold energy above the lowest energy one-magnon states, freezing out any potential decay channels for the lowest energy modes. In contrast, the experimentally observed spectrum appears to be broadened down to the lowest accessible energies suggesting that any gap is relatively small, if present at all. In fact, even in the extreme case where the spectrum is gapless, the decay matrix elements and density of decay states may be such that one-magnon intensity remains relatively sharp despite the presence of kinematically allowed one-magnon decay down to zero energy.

While the above remarks about the role of magnon interactions are largely model-independent, it is desirable to have a concrete calculation of their effects in Yb2Ti2O7 to establish the quantitative extent to which they capture the experimental data coming down from high fields. Speculations about the origin and nature of the observed continuum are moot in the absence of such a calculation so, in this paper, we directly address this issue. That this question can be meaningfully addressed in a quantitative manner rests on detailed fits of the sharp one-magnon dispersion measured in a regime of high fields where interaction effects can be neglected and LSWT is adequate to describe the physics Ross et al. (2011b); Thompson et al. (2017).

Due to the strong spin-orbit coupling of the $4f^{13}$ states of Yb3+, the symmetry-allowed exchanges can be highly anisotropic, as is indeed found in the parameters that have been determined empirically. We base our analysis on the parameter set of Ref. Thompson et al. (2017), determined by comparision to the spectrum in high magnetic fields. These exchange parameters place Yb2Ti2O7 close to the (semi-classical) phase boundary between the splayed ferromagnet (SFM) of Yb2Ti2O7 and an antiferromagnetic $\Gamma_{5}$ phase. It has been suggested that proximity to this phase boundary underlies the exotic properties of Yb2Ti2O7 Jaubert et al. (2015); Robert et al. (2015). This picture is consistent with the fact that hydrostatic pressure tunes samples that show no signs of long-range order into the SFM state Kermarrec et al. (2017).

The paper is organized as follows. First, in Sec. II, we describe the anisotropic exchange model for Yb2Ti2O7 considered in the remainder of the paper (Sec. II.1) and then outline our approach to NLSWT (Sec. II.2).

In Sec. III, we consider NLSWT corrections to static properties, in particular the SFM-$\Gamma_{5}$ phase transition and critical magnetic fields. In Sec. III.1, we study the zero field phase boundary using NLSWT, finding that the SFM phase becomes unstable near the boundary with the $\Gamma_{5}$ phase. We show that the parameters of Ref. Thompson et al. (2017), which are close to the classical phase boundary, sit almost exactly on top of this NLSWT instability line. This is our first main result: that the effects of magnon interactions drive Yb2Ti2O7 even closer to the SFM-$\Gamma_{5}$ phase boundary. In Sec. III.2, we study the phase transitions in applied magnetic fields along $[1\bar{1}0]$ and $[111]$ directions. Taking into account corrections from NLSWT, we find the critical fields depart from their classical values, which crucially bring them closer to what has been observed experimentally. This is our second main result.

We now turn to the key question: how does LSWT breakdown as the fields are lowered? These results are presented in Sec. IV, where we find that the regime of (semi-)quantitative reliability of spin-wave theory is significantly enhanced by the inclusion of magnon interactions. This is supported by a detailed study of the dynamical structure factor computed in NLSWT for three high symmetry field directions $[001]$ (Sec. IV.1), $[1\bar{1}0]$ (Sec. IV.2) and $[111]$ (Sec. IV.3). Next, in Sec. IV.4, we make a direct comparison with the existing experimental data of Thompson et al. (2017) for $[001]$ magnetic fields. For fields larger than $\sim 3\ {\rm T}$, where LSWT has been applied successfully, we confirm that departures from LSWT are small. For lower fields, we find that NLSWT captures well the renormalizations of the one-magnon modes down to fields of about $|\bm{B}|\sim 0.5\ {\rm T}$. For $[1\bar{1}0]$ fields we confirm the validity of LSWT for $|\bm{B}|\gtrsim 2\ {\rm T}$, as had been assumed in previous works Ross et al. (2011b). For $[111]$ fields, we find significant departures from LSWT for fields as high as $|\bm{B}|\sim 3\ {\rm T}$. In Sec. IV.5, we also provide a comparison to the high-resolution powder-averaged zero field inelastic data of Peçanha-Antonio et al. (2017), finding rough qualitative agreement.

In Sec. V, we summarize and discuss our results and their implications for Yb2Ti2O7. In Sec. VI, we discuss possible origins of the anomalous zero-field behavior, including intrinsic and extrinsic scenarios, as well as the importance of proximity to the SFM-$\Gamma_{5}$ phase boundary.

II Model and methods

We first describe the appropriate model for the magnetic ions in Yb2Ti2O7 and then outline our approach to NLSWT. Our treatment of magnon interactions is described briefly and we refer the reader to standard treatments for further details Chubukov et al. (1994); Chernyshev and Zhitomirsky (2009); Zhitomirsky and Chernyshev (2013).

II.1 Anisotropic exchange model

We consider the general nearest-neighbor model for effective spin-$1/2$ Kramers doublets on the pyrochlore lattice Curnoe (2008); McClarty et al. (2009); Onoda (2011); Onoda and Tanaka (2011); Ross et al. (2011b); Rau and Gingras (2018). One can restrict to the crystal field ground doublet in Yb2Ti2O7, as the gap to the next crystal field doublet Gaudet et al. (2015) is much larger than the exchange scale and the temperature. The appropriate effective model is then

[TABLE]

where the $\gamma_{rr^{\prime}}$ and $\zeta_{rr^{\prime}}$ are bond dependent phases and

[TABLE]

The relevant phases factors $\gamma_{ij}$ and $\zeta_{ij}$ for the three nearest-neighbor bond types (labelled $x$, $y$ and $z$, shown in Fig. 1) are

[TABLE]

where $\omega=e^{2\pi i/3}$ and the local axes, $(\bm{\hat{x}}_{r},\bm{\hat{y}}_{r},\bm{\hat{z}}_{r})$, follow the convention of Savary et al. (2012) (see App. A).

We primarily consider the parameters obtained in Thompson et al. (2017)

[TABLE]

where the Ising coupling is small, and $J_{z\pm}$ is the dominant exchange. These were obtained through fits to the magnetic excitations in high-fields in both $[1\bar{1}0]$ and $[001]$ directions Thompson et al. (2017) and are, further, consistent with the terahertz measurements of Ref. Pan et al. (2014). Earlier parametrizations, obtained using inelastic neutron scattering data at high field using a single field orientation, Ross et al. (2011b) and Robert et al. (2015), also lie in the classical ferromagnetic phase corresponding to the experimentally determined ground state. These differ, however, in detail, with the parameters of Ross et al. (2011b) being somewhat further from the classical phase boundary between the SFM and the $\Gamma_{5}$ states, than those of Thompson et al. (2017) or Robert et al. (2015).

II.2 Spin wave theory

To make our discussion somewhat general, we write the nearest-neighbor exchange model [Eq. (II.1)] in the form

[TABLE]

where $\bm{S}_{\bm{r}\alpha}$ is a spin-$S$ operator located in unit cell $\bm{r}$ with sublattice index $\alpha$. The spin operators can expressed in terms of Holstein-Primakoff bosons Holstein and Primakoff (1940); Kittel (1963); Auerbach (1998) as

[TABLE]

where $n_{\bm{r}\alpha}\equiv{a}^{\dagger}_{\bm{r}\alpha}{a}_{\bm{r}\alpha}$ and vectors $\bm{\hat{e}}_{\alpha,\pm}$, $\bm{\hat{e}}_{\alpha,0}$ define a local frame of reference with $\bm{\hat{e}}_{\alpha,0}$ being the classical ordering direction. It is useful to write the exchange matrix in this frame

[TABLE]

where $\mu,\mu^{\prime}=0,\pm$.

Expanding in powers of $1/S$ then yields a semi-classical expansion about the ordered state defined by $\bm{\hat{e}}_{\alpha,0}$. At $O(1/S^{2})$, this can be written as a sum of the classical energy, the usual two-magnon terms, as well as three- and four-magnon interactions. We write

[TABLE]

where the classical energy per site, $\epsilon_{\rm cl}$, is defined as

[TABLE]

and we have defined the individual pieces in symmetrized form as

[TABLE]

where $N_{s}$ is the number of sublattices, $N$ is the total number of sites and we have defined the Fourier transforms of the bosons as $a_{\bm{k}\alpha}\equiv N_{c}^{-1/2}\sum_{\bm{r}}e^{-i\bm{k}\cdot\bm{r}}a_{\bm{r}\alpha}$ where $N=N_{c}N_{s}$. Similarly, we have defined the Fourier transforms of the local exchange matrices as $\mathcal{J}^{\mu\mu^{\prime}}_{\bm{k},\alpha\alpha^{\prime}}\equiv\sum_{\bm{\delta}}\mathcal{J}^{\mu\mu^{\prime}}_{\bm{\delta},\alpha\alpha^{\prime}}e^{i\bm{k}\cdot\bm{\delta}}$. In terms of these local exchange matrices [Eq. (7)] one can write

[TABLE]

where the four-magnon vertices have been left unsymmetrized for brevity. As $1/S\rightarrow 0$, the interactions encoded in $H_{3}$ and $H_{4}$ can be included perturbatively as they are $O(S^{-1/2})$ and $O(S^{-1})$ with respect to the LSWT parts.

Mainly, we are interested in the dynamical structure factor, defined as

[TABLE]

where the real-time, Fourier transformed magnetic moment operators are defined $\bm{\mu}_{\bm{q}}(t)\equiv\sum_{r}e^{i\bm{q}\cdot\bm{r}}e^{iHt}\bm{\mu}_{r}e^{-iHt}$. A more useful form includes a polarization factor

[TABLE]

as this can be directly related to the intensity measured in an inelastic neutron scattering experiment Lovesey (1984) via

[TABLE]

where $F(\bm{q})\equiv F(|\bm{q}|)$ is the atomic form factor of Yb3+ Brown (2006).

The dynamical structure factor at $O(1/S^{2})$ requires the computation of the magnon Green’s function as well as several higher-order correlation functions. It is useful to consider three distinct pieces: the transverse-transverse part, which involves only the magnon Green’s function, $\bm{{G}}^{R}(\bm{k},\omega)$, the transverse-longitudinal parts, which involve three-magnon correlation functions, and the longitudinal-longitudinal parts which involve four-magnon correlation functions Mourigal et al. (2013). While the transverse-transverse part has $O(1/S)$ contributions, the other two parts appear first at $O(1/S^{2})$. We note that the Green’s function also appears in the transverse-longitudinal part of the dynamical structure factor, while the longitudinal-longitudinal part involves only the free magnon Green’s function at leading order Mourigal et al. (2013). Typically, the transverse-longitudinal and longitudinal-longitudinal parts are small relative to the leading transverse-transverse contributions. The central ingredient is then the (retarded) magnon Green’s function Ripka et al. (1986)

[TABLE]

where $\bm{\sigma}_{3}$ is a block Pauli matrix, $\bm{M}_{\bm{k}}$ is the (usual) LSWT dispersion matrix

[TABLE]

where $\bm{A}_{\bm{k}}={\bm{A}}^{\dagger}_{\bm{k}}$ and $\bm{B}_{\bm{k}}={\bm{B}}^{\intercal}_{-\bm{k}}$ are defined in Eqs. (11a) and (11b).

The self-energy, $\bm{\Sigma}^{R}(\bm{k},\omega)$, is generated by interactions, first at $O(S^{0})$. As for the Green’s function, this contains both normal and anomalous parts and sublattice indices. The relevant diagrams at $O(S^{0})$ in the magnon interactions are given in Fig. 2. We identify two distinct types of contributions to the self-energy: static (frequency independent) and dynamic (frequency dependent). The static contributions arise from Hartree-Fock-like diagrams involving the four-magnon interactions as well as tadpole-like diagrams arising from the three-magnon interaction. The dynamic contributions arise purely from the three-magnon terms. In addition to renormalizing the one-magnon spectrum they are also responsible for magnon decay Chubukov et al. (1994); Chernyshev and Zhitomirsky (2009), endowing the one-magnon states with finite lifetimes when kinematically allowed.

In practice, we first find the classical ground state numerically and compute the LSWT energies and eigenvectors on a grid of wave-vectors corresponding to a finite system of $N=4L^{3}$ sites with periodic boundary conditions. For static quantities, typically we find $L=8$ is sufficient to reach convergence, given the (generically) gapped spectrum (Sec. III). For dynamical quantities, we use $L=32$ when considering individual wave-vectors (Sec. IV.1, IV.2 and IV.3), and $L=16$ or $L=24$ when integrating over some range in momentum space (Sec. IV.4 and IV.5).

We then use these quantities to compute the required integrals in the self-energy diagrams numerically, working exclusively at zero temperature. These matrices then go into the Green’s function [Eq. (15)] along with a finite broadening $0^{+}\rightarrow\delta\sim O(10^{-3})$. The matrix inversion can then be carried out numerically and the Green’s function recovered; this effectively is a full solution of Dyson’s equation and avoids some of the issues discussed in Zhitomirsky and Chernyshev (2013) when trying to solve this equation perturbatively. This then yields the transverse-transverse part of the dynamical structure factor. Similar strategies are used for the transverse-longitudinal and longitudinal-longitudinal parts of the spin structure factor, with some (ad-hoc) resummations carried out in the $O(1/S^{2})$ corrections to the transverse-transverse and transverse-longitudinal parts to avoid introducing additional spurious poles, as described (for example) in  Mourigal et al. (2013).

Throughout we fix $S=1/2$, as appropriate for Yb2Ti2O7, extrapolating from the asymptotic limit $S\rightarrow\infty$ where spin-wave theory is controlled. Due to this choice our spin-wave expansion can break down, suggesting instabilities that appear only at small $S$, as we will see explicitly in Sec. III and Sec. III.2. This typically entails one of the (renormalized) excitation branches going soft, then continuing to negative (unphysical) energies. We also do not attempt to introduce self-consistency in the self-energy, be it in the static or dynamic parts; we will discuss some aspects of this in Sec. V. The results presented below, referred to as “NLSWT”, are computed following the above framework.

III Competing phases and renormalized phase boundaries

We begin with the effects of magnon interactions on the stability of the SFM state itself. Just as the softening of modes in LSWT suggest a classical instability, corrections to these phase boundaries can be obtained from NLSWT by considering where the (renormalized) spectrum goes soft. We explore these corrections in two contexts: in the location of the boundary between the SFM and $\Gamma_{5}$ phases in the anisotropic exchange model [Eq. (II.1)], and for the critical fields that separate the low-field SFM and the high-field polarized phase for two high symmetry field directions.

III.1 Splayed ferromagnet - $\Gamma_{5}$ phase boundary

We first consider the zero temperature phase diagram of the model [Eq. (II.1)] in the vicinity of the empirically determined exchange parameters of Thompson et al. (2017). These are more usefully written in the dual global frame (see App. A), in terms of the (dual) Heisenberg ($J$), Kitaev ($K$), symmetric off-diagonal ($\Gamma$) and Dzyaloshinskii-Moriya ($D$) exchanges Rau and Gingras (2018). In this language the parameters of Thompson et al. (2017) are

[TABLE]

For $J>0$, $D<0$ and $|K|,|\Gamma|\ll J,|D|$, as is relevant here, the phase boundary between the SFM and $\Gamma_{5}$ phase is $K+\Gamma\approx 0$ Rau and Gingras (2018, 2019). In Fig. 3 we show a section of the $(K+\Gamma)/J$-$D/J$ plane containing the parameters of Thompson et al. (2017). The classical phase diagram contains two phases - the SFM corresponding to the ground state of Yb2Ti2O7 and the $\Gamma_{5}$ phase that hosts an accidentally degenerate $U(1)$ manifold of states. Quantum fluctuations lift the $\Gamma_{5}$ degeneracy Jaubert et al. (2015); Yan et al. (2017) selecting the discrete set of $\psi_{2}$ or $\psi_{3}$ states (not shown) Poole et al. (2007). The empirically determined parameters for Yb2Ti2O7 lie close to this phase boundary as has been remarked upon previously Jaubert et al. (2015); Thompson et al. (2017); Yan et al. (2017).

Carrying out NLSWT about the classically stable SFM state shows that the lowest spin wave mode goes soft before reaching the classical phase boundary. This line of instability is shown in Fig. 3, reducing the range over which the SFM state appears in the $(K+\Gamma)/J$ direction. In the region between the line of instability and the classical phase boundary the state is a priori unknown. However, given the unstable mode appears with intensity at the $[111]$ wave-vector associated with the $\Gamma_{5}$ manifold, it is plausible that the phase that lies within this region is one of the order-by-quantum-disorder selected states $\psi_{2}$ or $\psi_{3}$. We stress that the line of instability potentially over-estimates the range of the SFM state, with a first order transition possible already at larger $(K+\Gamma)/J$. In addition, the effects of quantum fluctuations may introduce a finite gap to the soft mode, even at the phase boundary Rau et al. (2018).

We note that the parameters of Thompson et al. (2017) lie very close to the line of instability, just on the unstable (nominally $\Gamma_{5}$) side. This proximity reaffirms some of the arguments of Jaubert et al. (2015), where it was argued that the parameters of Ross et al. (2011b) were pushed closer to the SFM-$\Gamma_{5}$ phase boundary by quantum fluctuations. The calculations presented here suggest that this remains true for the parameters of Thompson et al. (2017); the proximity to this phase boundary likely plays an important role in understanding the physics of Yb2Ti2O7.

III.2 Critical magnetic fields

We now consider the behavior of the phase boundaries as a function of magnetic field. First, let us note that the instability of the parameters of Thompson et al. (2017) that appears at zero-field can be resolved by applying a very small $[001]$ magnetic field of the order of $\sim 10^{-3}\ {\rm T}$. As shown in Fig. 4, this opens a small gap for $|\bm{B}|\gtrsim 3\ {\rm mT}$ with intensity at the $[111]$ wave-vector. This can be contrasted with LSWT, where there is a finite $\sim 0.2\ {\rm meV}$ gap down to zero field. We thus see that, at small $[001]$ fields, there is a large renormalization of the $[111]$ mode, pushing it to low energy. The strong renormalization of this mode has mostly disappeared by $|\bm{B}|\sim 0.2\ {\rm T}$, as shown in Fig. 4. The smallness of the required stabilization field reflects the extreme proximity of these parameters to the line of instability in the phase diagram discussed in Sec. III.1, and shown in Fig. 3. We note that this zero-field instability is not present in the parameters of Ross et al. (2011b), due to its distance from the boundary, but is present for the parameters of Robert et al. (2015), which are (in some sense) closer to the classical boundary than those of Thompson et al. (2017) (and thus are in the region where NLSWT breaks down).

In contrast to the $[001]$ field direction, there are field-induced phase transitions out of the SFM state for both $[1\bar{1}0]$ and $[111]$ field directions, due to their remnant two-fold symmetry. To analyze what effects magnon interactions have on these phase boundaries, we track the evolution of the lowest lying mode at the $[000]$ wave-vector as a function of field. We restrict ourselves to the parameters of Thompson et al. (2017) throughout.

The result for $[1\bar{1}0]$ field is shown in Fig. 5(a). We see that coming down from the high-field phase the LSWT gap closes at $|\bm{B}|\sim 0.84\ {\rm T}$, well above the $|\bm{B}|\sim 0.5\pm 0.125\ {\rm T}$ 111The quoted error bars are our estimate due to the experiments being carried out for only a few fields near the critical field, namely at $|\bm{B}|=0.25\ {\rm T}$, $0.5\ {\rm T}$ and $0.75\ {\rm T}$ Ross et al. (2009). determined experimentally Ross et al. (2009) . The NLSWT result breaks down artificially above the classical field, due to the softening the bare LSWT spectrum. However from Fig. 5(a) it is clear that the gap at $[000]$ is strongly enhanced, suggesting that the critical field is lowered by the effects of magnon interactions, as is found experimentally. Given that NLSWT is limited by the classical phase boundaries, it is not directly possible to obtain a lowering of the critical field in this fashion. To resolve this, we have also looked at this gap coming out of the ordered phase. We find that in NLSWT the SFM state has an instability at $[000]$ which goes soft at $|\bm{B}|\sim 0.61\ {\rm T}$ [see Fig. 5(a)], somewhat above the experimental value but not quantitatively inconsistent.

The result for $[111]$ field is shown in Fig. 5(b). The classical transition field in this case, $|\bm{B}|\sim 0.46\ {\rm T}$, is lower than the $|\bm{B}|\sim 0.65\ {\rm T}$ that has been recently observed experimentally Scheie et al. (2017). Both theoretically and experimentally this is expected to be a first-order transition, with the gap remaining open across the transition. In contrast to the $[1\bar{1}0]$ case it should thus be possible to obtain an increase in the transition field from NLSWT in the polarized phase. One finds that the gap at $[000]$ is strongly reduced by interaction effects, finally going soft very close to the experimental value of $|\bm{B}|\sim 0.65\ {\rm T}$. Given the first order nature of the transition seen experimentally, we expect the precise matching of these results to be somewhat accidental, with the true (theoretical) transition being first order at a field above the NLSWT instability.

IV Dynamical Structure Factor

We now explore more broadly the interaction corrections from NLSWT to the dynamical structure factor. We consider both the parameter sets of Ross et al. (2011b) and Thompson et al. (2017) as well as fields in the $[001]$ (Sec. IV.1), $[1\bar{1}0]$ (Sec. IV.2) and $[111]$ (Sec. IV.3) directions. To obtain a global picture of the spectrum, we compute the dynamical structure factor along paths in momentum space between high symmetry points, that (mostly) lie in the plane perpendicular to the applied magnetic field. The results are shown in Fig. 6 ($[001]$ field), Fig. 7 ($[1\bar{1}0]$ field) and Fig. 8 ($[111]$ field). For each we show the full dynamical structure factor, including the transverse-transverse, transverse-longitudinal and longitudinal-longitudinal components (as discussed in Sec. II.2). This result is convolved with a narrow Gaussian to improve visibility of any sharp modes. In addition to the structure factor in NLSWT, we also include the bare LSWT one-magnon spectrum and the edges of the bare two-magnon continuum.

After this overview, we consider two direct comparisons to experimental data, namely the single crystal data of Thompson et al. (2017) (Sec. IV.4) and the powder-averaged data of Peçanha-Antonio et al. (2017) (Sec. IV.5). For each case the integration over momentum space performed experimentally is reproduced theoretically with the form factor included [Eq. (14)].

IV.1 [001] magnetic field

First, let us revisit the dependence of the excitations on a $[001]$ field, as was used by Thompson et al. (2017) to determine the exchanges and $g$-factors. A sequence of fields, from $|\bm{B}|=5\ {\rm T}$ down to $|\bm{B}|=0^{+}\ {\rm T}$ (to select a single domain) is shown in Fig. 6 using the parameters of Thompson et al. (2017). Examining the LSWT and NLSWT results at $|\bm{B}|=5\ {\rm T}$ (see Fig. 6(a)), we see that the two-magnon continuum is well-separated from the one-magnon states, and corrections due to magnon interactions are small. We conclude that LSWT is sufficient to describe the spectrum for $|\bm{B}|\gtrsim 5\ {\rm T}$, and thus the fitting procedure used in Thompson et al. (2017) is sensible, up to the level of precision implied in Fig. 6(a).

As we lower the field, interaction effects become important, with significant differences between the LSWT and NLSWT spectrum becoming apparent for $|\bm{B}|\lesssim 3\ {\rm T}$. These are the fields at which the two-magnon continuum begins to overlap with the one-magnon bands [see Fig. 6(c)]. Roughly speaking, one may rationalize the renormalizations as a kind of level repulsion between the one-magnon spectrum and the encroaching two-magnon continuum. When the one- and two-magnon states overlap, spontaneous magnon decay is kinematically allowed. Such decay is most apparent for fields $|\bm{B}|\lesssim 2\ {\rm T}$ near the upper bands at $[2\bar{2}0]$, as shown in Figs. 6(d-f). In such regions, the one-magnon excitations acquire a true linewidth, as opposed to the artificial width given to the (sharp) magnon modes outside these regions. As we lower the field further, this mode becomes more and more damped, as the two-magnon density of states increases at its energy.

Below $|\bm{B}|\lesssim 1\ {\rm T}$ the spectrum begins to stabilize with only small, quantitative changes appearing, as shown in Fig. 6(e,f). In contrast to experiments, this calculation exhibits sharp magnon modes across most of momentum space, within the energy range $0.2\ {\rm meV}\lesssim\omega\lesssim 0.7\ {\rm meV}$, save for a few regions where decay is kinematically allowed. We further note that the intense magnon mode at $0.2\ {\rm meV}$ near $[000]$ remains present (at roughly the same energy as in LSWT), but is significantly flattened. Indeed, the NLSWT corrections remove nearly all of the dispersion along the $[h00]$ direction, as can be seen in Fig. 6(e). Similar magnon band flattening has been observed in interaction corrections to the spin one-half triangular lattice antiferromagnet Zheng et al. (2006); Starykh et al. (2006). In addition to these sharp one-magnon modes, the intensity of the two-magnon continuum is also enhanced at low fields, appearing as broad intensity in the range $0.5\ {\rm meV}\lesssim\omega\lesssim 1.0\ {\rm meV}$.

Finally, as discussed in Sec. III.2, at strictly zero-field the SFM state is unstable, with a mode going slightly soft at $[111]$. Remnants of this mode can be seen near $[\bar{2}00]$ and $[020]$ in Fig. 6(f). By application of a small field ($|\bm{B}|\sim 3\ {\rm mT}$) this mode can be lifted to finite energy. The qualitative features of the spectrum discussed above are not strongly affected by the presence or absence of this mode, given its low density of states. However, it would appear as an intense, nearly gapped mode near $[111]$, signalling the proximity to the SFM-$\Gamma_{5}$ phase boundary (see Sec. V for further discussion).

IV.2 $[1\bar{1}0]$ magnetic field

Next we consider the dynamical structure factor for an applied $[1\bar{1}0]$ field, as shown in Fig. 7. We focus on the two experimentally studied fields $|\bm{B}|=2\ {\rm T}$ and $|\bm{B}|=5\ {\rm T}$, as first studied in Ref. Ross et al. (2011b), both above the critical field in the polarized phase, using the parameters of Thompson et al. (2017).

As in the case of $[001]$, at $|\bm{B}|=5\ {\rm T}$ the NLSWT calculation matches LSWT almost perfectly, further justifying the use of LSWT in fitting the model at these fields [Fig. 7(a)]. The applicability of LSWT at the lower field, $|\bm{B}|=2\ {\rm T}$, is less clear on naïve grounds. From simple kinematic considerations one may worry about interaction effects becoming important, as the two-magnon continuum has begun to overlap with the one-magnon states, as shown in Figs. 7(b). However, while there is some decay and some renormalization, its effects are rather small relative to (say) typical experimental resolution, and so for $[1\bar{1}0]$ even LSWT is likely sensible at $|\bm{B}|\gtrsim 2\ {\rm T}$. We have confirmed that this is somewhat generic, also holding for the parameters of Ross et al. (2011b) (not shown).

IV.3 $[111]$ magnetic field

We now turn to results for the third high symmetry field direction $[111]$, which has been of interest recently Scheie et al. (2017). Here we find considerable departures from LSWT even at $|\bm{B}|=3\ {\rm T}$ [see Fig. 8(a)], well above the transition field of $|\bm{B}|\sim 0.65\ {\rm T}$ Scheie et al. (2017). Even at this high field several of the one-magnon bands strongly intersect with the two-magnon continuum. The highest energy mode, originating from the pinned moment of the single sublattice distinguished by the field direction, lies entirely within the two-magnon continuum. This induces renormalization of the energies as well as spontaneous magnon decay, for example near [$2\bar{2}0]$ as shown in Fig. 8(a).

Going down to $|\bm{B}|=1.5\ {\rm T}$, still significantly above the transition field of $|\bm{B}|\sim 0.65\ {\rm T}$ Scheie et al. (2017), there are drastic corrections to the LSWT spectrum. While lowest mode remains sharp, it receives large corrections, especially where it meets the second and third modes of the spectrum. These two modes are nearly completely destroyed by spontaneous magnon decay; over most of the wave-vectors shown in Fig. 8(b) there is either significant decay or significant renormalization, or both. The $[111]$ polarized phase thus appears to be an excellent setting to experimentally observe and study spontaneous magnon decay. We note that the highest energy band is relatively sharp and only weakly renormalized by magnon interactions. Finally, we have also confirmed that the features discussed above are also present in the parameters of Ross et al. (2011b) (not shown), and thus appear to be relatively generic for $[111]$ fields. However, given the small one-magnon gap at $1.5\ {\rm T}$, and thus the low-lying two-magnon continuum, details of the renormalization and induced decay do depend somewhat on the precise exchanges considered.

IV.4 Direct comparison to experiment in [001] field

We now consider a more direct comparison between the LSWT and NLSWT dynamical structure factor and the detailed experimental data in $[001]$ magnetic field, as studied in Ref. Thompson et al. (2017). This field direction selects a single domain at low fields and adiabatically connects the low field ferromagnet to the polarized state. We focus on three cuts in momentum space: one along $[h00]$ (Fig. 9), one along $[\bar{\zeta}{\zeta}0]$ (Fig. 10) and one along $[-1+\xi,1+\xi,0]$ (Fig. 11). For each cut we integrate over a region out-of-plane and perpendicular to the cut line of width $[-0.2,0.2]$ r.l.u to allow direct comparison to experiment Thompson et al. (2017). We compute the full neutron scattering intensity, including the Yb3+ atomic form factor, as given in Eq. (14). Throughout we use the parameters of Ref. Thompson et al. (2017) and fix the overall intensity scale globally, for all fields and cuts, to best match the $[h00]$ cut at $|\bm{B}|=5\ {\rm T}$.

To frame our discussion, let us highlight three features of the data and calculations; first, we have the one-magnon-like features that are sharp in energy and wave-vector with high-intensity. Second, we have the extent of broadening of the magnon bands. Third, we have the intensity and location of the (primarily) two-magnon continuum. We note that, as the two-magnon continuum is not determined self-consistently in our NLSWT, it will be systematically higher in energy than what may be expected from the one-magnon energies, which are renormalized downward.

As in the previous subsection, we find that the results at $|\bm{B}|=5\ {\rm T}$ match the experimental data very well, with NLSWT offering no significant improvement over LSWT [see Figs. 9(a), 10(a), 11(a)]. We note that the agreement of one-magnon energies from NLSWT are slightly worse than those from LSWT [see, e.g., Fig. 10(a)]; this is to be expected given that these parameters were determined using LSWT and there are small interaction corrections even at these high fields.

At $|\bm{B}|=1.5\ {\rm T}$, the field is sufficiently low that the one- and two-magnon states now overlap. One finds qualitative agreement between NLSWT and LSWT for primarily one-magnon features, but with visible quantitative differences in the intensity of the two-magnon continuum, and in some broadening [e.g. near $[\bar{1}00]$ Fig. 9(b)]. Across each of the cuts these differences improve the agreement of NLSWT with experiment relative to LSWT. NLSWT also includes some intensity from the two-magnon continuum, roughly in the same location as some of the broad, diffuse scattering seen experimentally [shown in Figs. 9(b), 10(b), 11(b)], though it is somewhat weaker than in the theoretical calculation.

At lower fields neither LSWT nor NLSWT give a full account of the experimental data. The experimental result includes a very diffuse background of scattering over a wide range in energies that is neither reproduced in NLSWT by spontaneous magnon decay nor by the direct intensity of the two-magnon continuum. However, many of the one-magnon features of NLSWT match well with the sharper features of the experimental data. For example, the nearly flat mode along $[h00]$ at $|\bm{B}|=0.21\ {\rm T}$ [see Fig. 9(c)] that is observed experimentally is also seen in NLSWT, but not in LSWT where it has significantly more dispersion. We find similar rough agreement also for the location of the gross features of the $[\bar{\zeta}{\zeta}0]$ and $[-1+\xi,1+\xi,0]$ cuts [see Figs. 10(c), 11(c)]. As in the $1.5\ {\rm T}$ data the intensity of the two-magnon continuum is also markedly lower in NLSWT than in the experimental data. In the zero-field result, this trend is enhanced, with the experimental data becoming broader and more diffuse, but with the gross intensity features being in roughly the same locations as in NLSWT [see Figs. 9(d), 10(d), 11(d)].

IV.5 Powder-averaged inelastic neutron scattering

Given the sample dependence present in crystals of Yb2Ti2O7 grown under different conditions or by different groups, it is worthwhile to make some comparisons to other data sets. Here we consider the inelastic data on the powder of Peçanha-Antonio et al. (2017), as it has higher energy resolution and shows more structure than (similar) powder data obtained by other groups. Other powder samples, for example as presented in the review of Hallas et al. (2018), show only a broad continuum, similar to what is seen in the single-crystal of Thompson et al. (2017). Since we have compared to this data at length in Sec. IV.4, we do not make direct comparisons to these other powder samples here.

For a zero-field cooled sample, powder-averaging can be carried out straightforwardly, by integrating over wave-vectors with a given magnitude. We present the powder averaged intensity [Eq. (14)] for the parameters of Ref. Thompson et al. (2017) in Fig. 12. We have added a small ($10\ {\rm mT}$) $[001]$ field to render the SFM state stable. For comparison we also include the data of Peçanha-Antonio et al. (2017) presented similarly. We see that several of the gross features are reproduced: a gap and band of intensity near $\sim 0.2\ {\rm meV}$ which is maximal at small $|\bm{Q}|$ and a broad continuum of intensity extending from $\sim 0.3\ {\rm meV}$ out to $\sim 0.8\ {\rm meV}$ with maximal intensity near $\sim 1\AA^{-1}$. Several features are absent however, most notably the high intensity near $1.1\AA^{-1}$ and the sub-gap features visible in the range $0.5\AA^{-1}\lesssim|\bm{Q}|\lesssim 1.1\AA^{-1}$. We do note that this intensity near $1.1\AA^{-1}$ coincides with location of the nearly gapless mode at $[111]$ that is present in our calculation. However the intensity of this mode does not match; as the minimum only occurs in a small region of momentum space, only a weak intensity is visible extending down at $1.1\AA^{-1}$ in the NLSWT calculation. As in Sec. IV.4, the intensity of the higher energy two-magnon continuum is underestimated by the NLSWT.

V Discussion

We have presented a detailed study of NLSWT as applied to Yb2Ti2O7. In the fully polarized phase where fluctuations are small and multi-magnon continua are well-separated from the one-magnon states, LSWT is expected to provide a quantitative description of the one-magnon dispersions. This observation motivated fits of LSWT to inelastic neutron scattering data at $5\ {\rm T}$ and above where the magnetic field energy is comfortably greater than the exchange scale. Our NLSWT calculations confirm that the effects of magnon interactions are weak at these fields providing a retrospective justification for the LSWT fits.

The corrections from NLSWT are expected to be more significant at lower fields, and to be quantitative within some window of intermediate fields. Indeed, we find that the critical fields predicted by LSWT are significantly renormalized by magnon interactions in the direction expected on the basis of experiments (Sec. III.2). By comparing with inelastic neutron scattering data (Sec. IV), we have found empirically that one-magnon energy renormalizations and spontaneous magnon decay within the higher energy magnon bands starts to become significant for $[001]$ and $[1\bar{1}0]$ fields $\lesssim 2\ {\rm T}$ and for $[111]$ fields $\lesssim 3\ {\rm T}$. Comparison to Ref. Thompson et al. (2017), which contains energy-momentum slices for various $[001]$ fields, suggests that NLSWT is nearly quantitative at $\sim 1.5\ {\rm T}$, and captures some aspect of the experimental intensity reasonably well down to even $\sim 0.2\ {\rm T}$ (see Figs. 9, 10 and 11), but misses the extensive broadening. This is a central result of this paper.

We turn now to the low field data. Inelastic neutron scattering data at zero field can be found as powder averages in Refs. Gaudet et al. (2016); Peçanha-Antonio et al. (2017) and in single crystal samples in Refs. Robert et al. (2015); Thompson et al. (2017). In Refs. Robert et al. (2015); Gaudet et al. (2016); Thompson et al. (2017) the scattering intensity appears to have no sharp features within experimental resolution, though clear intensity modulations are seen within the apparent continuum; see Figs. 9(d), 10(d) and 11(d) for reproductions of the Ref. Thompson et al. (2017) zero field data. The powder data of Peçanha-Antonio et al. (2017) is compared with the ($10\ {\rm mT}$) NLSWT calculation using the parameters of Thompson et al. (2017). While continuum scattering is visible in that data above $\sim 0.3\ {\rm meV}$, the low energy sharp features present in the calculation, such as the flat band of intensity at about $\sim 0.2\ {\rm meV}$, are present in the data.

In the calculation, there is an additional soft mode with intensity at $[111]$. We note that there is a corresponding feature in the data of Peçanha-Antonio et al. (2017), albeit with a higher intensity than predicted by NLSWT. While some evidence has been presented Chang et al. (2016) for a $\sim 40\ \mu{\rm eV}$ soft mode at $[111]$ in a single crystal via inelastic neutron scattering at low fields, this has not been explored in detail in the literature. That this gap has not been seen in other crystals may simply be a consequence of the smallness of the gap relative to the energy resolution of these other experiments.

On the theoretical side, the presence of such a soft mode is expected near the phase boundary between the SFM ground state of Yb2Ti2O7 and the antiferromagnetic $\Gamma_{5}$ manifold Jaubert et al. (2015); Yan et al. (2017), as has been argued to be relevant for the exchange parameters of Refs. Ross et al. (2011b); Robert et al. (2015); Thompson et al. (2017). The phase boundary obtained classically is renormalized by quantum fluctuations Jaubert et al. (2015), and we find that the parameters of Ref. Thompson et al. (2017), noted as being close to the classical phase boundary, are even closer to the renormalized instability line. These observations suggest that Yb2Ti2O7 may be closer to the SFM-$\Gamma_{5}$ phase boundary than would be expected on the basis of a classical calculation from the empirically determined exchange parameters.

Furthermore, given the presence of a gapped mode at the zone center, one can make the model-independent observation that naïve kinematics already forbid broadening of this mode via magnon interactions. Leading order NLSWT, as we have presented here, has difficulty capturing this quantitatively, as the two-magnon continuum is not determined self-consistently and reflects the bare one-magnon kinematics. However, if a self-consistent treatment that resolves such issues were carried out, it would likely not be able to capture the decay region seen experimentally, given it still lies outside the (loose) bounds defined by even the experimental one-magnon features. These conclusions are unaffected by the presence or absence of the nearly soft mode at $[111]$. While such a soft mode can place much of the rest of the one-magnon spectrum in the two-magnon continuum, the density of states associated with this mode is too small to induce any significant linewidth. One potential resolution could be in a breakdown of the naïve kinematics described above through strong interactions within the two-magnon continuum Zhitomirsky and Chernyshev (2013); Verresen et al. .

VI Outlook

We have shown that NLSWT describes the single crystal inelastic neutron scattering data at intermediate strength $[001]$ magnetic fields, capturing some of the observed broadening of the one-magnon lineshapes and the renormalization of the one-magnon energies. At lower fields, the theory fails to account for the anomalously broadened magnetic scattering observed in various powder and single crystal samples. There are several possible scenarios for how this could be resolved, which we loosely categorize as intrinsic and extrinsic, i.e. based on a more detailed analysis of the present model, or involving ingredients beyond it.

Several possible intrinsic effects, beyond the magnon interactions we have considered here, could play a role in the low-field physics of Yb2Ti2O7. We have described ways in which we expect our calculations to fall short of capturing the full extent, at low fields, of the interaction corrections to LSWT. Further investigation of the renormalization of the two-magnon states and their self-consistent influence on the one-magnon branches may be useful in this regard. We stress, however, that for this kind of explanation to capture the zero field excitations in Yb2Ti2O7, the naïve kinematics outlined in Sec. V must break down to allow the broadening to persist to the lowest energies, as observed experimentally.

Also, it is possible that multi-magnon effects involving three or more magnons may be required to understand the anomalous broadening. In such a case it may make sense to describe these excitations, while still arising from the conventional ferromagnetic phase, as being related to some nearby exotic “parent” state, such as some kind of quantum spin liquid Ross et al. (2011b); Chang et al. (2012); Sanyal et al. ; Chern and Kim . Similar scenarios have been outlined to understand the high-energy excitations in the Kitaev magnet $\alpha$-RuCl3, with the Kitaev spin liquid serving as the parent state Banerjee et al. (2016). Given the proximity to the SFM-$\Gamma_{5}$ phase boundary, it is not implausible that quantum effects could stabilize a disordered intermediate phase near the boundary Sanyal et al. ; Chern and Kim that could serve as parent state.

While intrinsic effects are clearly important at low fields, the strong sample to sample variations that have been reported Ross et al. (2012); Gaudet et al. (2015); Mostaed et al. (2017); Shafieizadeh et al. (2018); Bowman et al. (2019) in Yb2Ti2O7, as well as the issues raised in Sec. V, suggests that extrinsic effects may also be important at low fields. This calls into question whether further refinements to spin wave theory are capable of explaining all of the puzzling aspects of the zero field data. This sample dependence, along with the sensitivity to applied hydrostatic Kermarrec et al. (2017) or chemical Hallas et al. (2018) pressure, further supports an extrinsic origin of some of the broadening exhibited at low fields. If, indeed, extrinsic effects contribute significantly to the scattering continuum, sufficiently clean crystals may exhibit inelastic spectra more in line with predictions of NLSWT. From this perspective, the relatively sharp features in the zero field inelastic powder data of Peçanha-Antonio et al. (2017) and the sharp low energy soft mode reported in Chang et al. (2016) might be interpreted as arising from the underlying intrinsic physics.

Extrinsic effects, such as structural disorder, can take many distinct forms. Several different types of structural disorder have been discussed for Yb2Ti2O7: these include point-defects, such as “stuffed” Yb3+ spins Ross et al. (2012); Gaudet et al. (2015); Bowman et al. (2019), oxygen vacancies Mostaed et al. (2017); Bowman et al. (2019), as well as internal strains and atomic displacements, such as those that might arise due to line-defects like dislocations Shafieizadeh et al. (2018) or through intrinsic frustration of the structural network Trump et al. (2018). For example, recent experiments have explored the effects of oxygen annealing and stuffed Yb3+ spins Bowman et al. (2019), finding that the presence or absence oxygen vacancies has strong effects on the magnetic correlations. It has also been reported Shafieizadeh et al. (2018) that some single crystals are liable to contain defects such as “superdislocations” or anti-phase boundaries that can induce further lattice distortions.

Any of these defects could provide a route to broadening the spectrum. However, given the (nominally) low-level of such defects reported in many samples Bowman et al. (2019); Arpino et al. (2017), one then must grapple with the question of how (nominally) small amounts of disorder could lead to the large observed broadening. Indeed, the sensitivity Arpino et al. (2017) to such extrinsic effects suggests that the physics of Yb2Ti2O7 may be quite delicate. One potential route to rendering the magnetic physics delicate, and thus extremely sensitive to disorder, is close proximity to the SFM-$\Gamma_{5}$ boundary Jaubert et al. (2015); Robert et al. (2015). One could then imagine that defects disturb the local physics enough to mix the SFM and $\Gamma_{5}$ states, strongly affecting the nature of the dynamical response.

As a concrete example, one can consider domain walls between the different SFM domains. Near the boundary between the SFM and $\Gamma_{5}$ phases, these domain walls in fact (locally) pass through the manifold of $\Gamma_{5}$ states, reflecting the topology of the piecewise connected, (nearly) degenerate $U(1)$ manifolds Canals et al. (2008); Chern ; Yan et al. (2017). The width of these domain walls diverges (classically) as one approaches the boundary, potentially providing a route to significant broadening even for relatively low levels of disorder. 222We note that such a possibility, with the $\Gamma_{5}$ states appearing in SFM domain walls, was recently independently noted by the Johns Hopkins group Scheie et al. (2019)

Crucially, as we showed in Sec. III.1, the effects of magnon interactions appear to move Yb2Ti2O7 even closer to the phase boundary, further enhancing any such effects. Such a scenario would also be largely compatible with the dynamics being simpler at high-fields: once the magnetic field energy lifting the near-degeneracy at low fields dominates the disorder energy scales, the phase would simply be a SFM with a small amount of disorder, and thus would likely exhibit a significantly smaller amount of broadening than when the SFM and $\Gamma_{5}$ states are mixed. With both intrinsic and extrinsic scenarios in contention, there remains a great deal to do to reach a quantitative understanding of the zero field spectrum in Yb2Ti2O7. Given the progress in understanding obtained starting from the high field limit, we are optimistic about the prospects for settling this question in Yb2Ti2O7 in the near future.

The importance of resolving this issue is highlighted by the the existence of puzzling continua in other magnetic materials. While in the foregoing discussion we have focused on the experimental results in Yb2Ti2O7, there are tantalizing results Hallas et al. (2018) that suggest that the physics seen in Yb2Ti2O7 also appears in some other, less explored, materials such as Yb2Sn2O7 and Yb2Ge2O7. These two compounds are isostructural to Yb2Ti2O7, with the substitution of the non-magnetic Ti for Sn or Ge expected mainly to act as chemical pressure Dun et al. (2014). One, Yb2Sn2O7, is known to have a SFM ground state Dun et al. (2013). while Yb2Ge2O7 Hallas et al. (2016) shows $\Gamma_{5}$ antiferromagnetic order. However, both show a similar broad continuum of excitations Hallas et al. (2016, 2018) at zero-field. While a detailed parametrization of the exchange in Yb2Sn2O7 or Yb2Ge2O7 has not yet been carried out, it is plausible that they also live close to the SFM-$\Gamma_{5}$ phase boundary Jaubert et al. (2015), like Yb2Ti2O7. In the same vein, it has been argued Dalmas de Réotier et al. (2017) that the rare-earth spinels AYb2X4 (A=Cd, Mg and X=S,Se) Lau et al. (2005); Higo et al. (2017); Dalmas de Réotier et al. (2017), which exhibit $\Gamma_{5}$ order Dalmas de Réotier et al. (2017), could share much of the physics. However, given their different crystal structure, extrinsic effects could manifest themselves differently and may shed some light on the physics of Yb2Ti2O7, Yb2Ge2O7 and Yb2Sn2O7.

Acknowledgements.

We thank R. Coldea and M. Gingras for helpful comments and collaborations on prior work. We also thank R. Coldea as well as V. Peçanha-Antonio and Y. Su for providing experimental data (Refs. Thompson et al. (2017) and Peçanha-Antonio et al. (2017) respectively). This work was in part supported by Deutsche Forschungsgemeinschaft (DFG) under grant SFB 1143 and through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490).

Appendix A Conventions

We choose the four local frames $(\bm{\hat{x}}_{i},\bm{\hat{y}}_{i},\bm{\hat{z}}_{i})$ for the pyrochlore lattice following Savary et al. (2012)

[TABLE]

where $\bm{\hat{y}}_{i}=\bm{\hat{z}}_{i}\times\bm{\hat{x}}_{i}$. The four basis sites of the tetrahedron are then along the $\bm{\hat{z}}_{i}$ directions, with $\bm{\hat{r}}_{i}=\bm{\hat{z}}_{i}$.

Exchange parameters have been reported in the main text in this local frame, via $J_{zz}$, $J_{\pm}$, $J_{\pm\pm}$ and $J_{z\pm}$ and in the dual global frame Rau and Gingras (2018) where the symmetry-allowed couplings on one pair of nearest neighbors (sublattices $1$ and $2$) can be written

[TABLE]

with transformation

[TABLE]

between the exchange couplings. Here $J$, $K$, $D$, $\Gamma$ are respectively couplings for (dual) Heisenberg, Kitaev Dzyaloshinskii-Moriya, and symmetric off-diagonal exchange couplings.