Inelastic neutron scattering data

We first examine the low-energy magnetic excitation spectra of YbCl₃ at 0.24 K. Figures 2a–c and 3a–e show the INS spectra as a function of energy, ℏω, and wave vector, Q, transfer. Figure 2a–f is plotted as the product of the intensity and the energy transfer to emphasize higher energy features in the spectrum. The spectra contain three distinct features: a component characteristic of conventional transverse spin waves (ℏω ≤ 0.6 meV), a continuum or multimagnon component, and a sharper component at higher energies (0.8 ≤ ℏω ≤ 1.2 meV). The spin-wave mode disperses throughout the (HK0) plane with no discernible dispersion along the (00L) direction (see Fig. 3e and Fig. 5 of the Supplementary material). The lack of an observable dispersion along the (00L) direction indicates that the interactions between the honeycomb planes are very weak. Interestingly, the spin-wave mode appears to be broader than the instrumental resolution. While the reason for this broadening is uncertain, it may arise due to domain formation as a result of bond disorder⁶⁵, interaction with the multimagnon continuum (potentially following a renormalization of the spin-wave spectrum⁵²), or even the inability to fully distinguish the overlapping continuum contribution. We also point out that the T = 3.7 K data in Fig. 3f illustrate a complete lack of well-formed magnetic excitations at higher temperatures. Finally, an important feature of the data is the lack of an appreciable spin gap (see Supplementary material for additional details). This observation suggests that the spins do not possess a significant uniaxial anisotropy, in agreement with the crystal field ground state with easy plane anisotropy determined in ref. ⁶¹.

Fig. 2

Inelastic neutron scattering data and linear spin-wave theory calculations.

a–c Inelastic neutron scattering data measured at T = 0.24 K using the CNCS instrument. d–f Linear spin-wave calculations including transverse and longitudinal channels for the ideal model (see Eq. (1)) with a fitted value of J = 0.421(5) meV along high symmetry directions in the (HK0) plane. Panels a–f are plotted as the product of intensity and energy transfer (ℏω). The Yb³⁺ magnetic form factor is included in the calculations. g Intensity as a function of ℏω through the Van Hove singularities at several wave vectors for T = 0.24 K (T = 5 K), marked by solid (open) points. Solid lines are Gaussian lineshapes with a sloping background. The horizontal bar represents the energy resolution at ℏω = 1.19 meV. The color of the data in (g) corresponds to the wave vector indicated by the colored arrow at the top of (a)–(c). Peak positions of the Gaussian lineshapes in (g) are shown as solid circles in (d)–(f). h–i Calculated (left) and measured (right) scattering intensity for ℏω = 0.9 meV (h) and 0.3 meV (i). The data and calculation in (h) have been scaled by a factor of 5 to be on the same intensity scale as (i). Gray lines indicate high symmetry directions of the Brillouin zone. Wave vector transfers are shown in the reciprocal space of the monoclinic lattice in reciprocal lattice units for (a)–(g) and projected into a hexagonal lattice in (h)–(i). Data in (a)–(c) and (g) have been averaged about the H = 0 and K  = 0 axes. Data in (h) and (i) have been averaged about the vertical axis.

Fig. 3

Inelastic neutron scattering data and calculations for YbCl₃.

a–e Measured at 0.24 K using the CNCS instrument. f Measured at 3.7 K. Black points are the locations of the absolute peak intensity at different wave vectors. The solid black line is the fitted spin-wave dispersion with a single interaction J = 0.421(5) meV. The dashed red line is the calculated upper boundary of the corresponding two-magnon continuum scattering. g–k Linear spin-wave calculations of the transverse spectrum for interactions J = 0.42(1) meV, $$ J^{\prime }=0.43\left(1\right)$$ meV, and J_c = 0 meV, as described in the text. Data in (a)–(f) have been averaged about the H = 0 and K = 0 axes.

The most unusual part of the spin excitation spectrum is the sharp feature toward the top of the broad continuum. While there is precedence for the observation of spinon and multimagnon continua in one-dimensional^66–70 and two-dimensional square and triangular lattice^55–60 quantum magnets, the observation of a sharp feature within such a continuum has, to the best of our knowledge, not yet been reported. This sharp multimagnon feature is explored further through constant wave vector scans in Fig. 2g. At energies where the sharp feature is well separated from the most intense parts of the continuum, the width is essentially limited to the calculated energy resolution of the instrument, FWHM = 0.04 meV at ℏω = 1.19 meV. As will be described in more detail below, this result provides an experimental demonstration of how a sharp feature can be generated within a continuum response in a quantum magnet. Note that the hexagonal symmetry of the spin excitations is manifest in both these higher energy features and the transverse spin-wave modes at lower energies, as shown in the right sides of Fig. 2h and i, respectively.

Linear spin-wave theory

To understand the physics begetting the novel spin excitation spectrum of YbCl₃, we consider a Heisenberg model on the honeycomb lattice with a single antiferromagnetic exchange interaction J between nearest-neighbor S =  1/2 spins

Comparison between data and model

The HLHM in Eq. (1) reproduces the experimental data for the transverse spin-wave mode, the broad continuum, and the sharp feature toward the top of the continuum. We first note that, due to the summation over the momentum k, the contribution from the two-magnon states in Eq. (5) results in a broad continuum of scattering. To determine the Heisenberg exchange J, we consider the transverse component of the data. Due to the large scattering intensity of the continuum, we compare the calculated dispersion to the overall maxima in the scattering intensity as a function of q (solid points in Fig. 3a–e). Comparing these values for points restricted to the (HK0) plane yields a Heisenberg exchange of J = 0.421(5) meV, shown as a solid line in Fig. 3a–e. Moreover, we directly compare the data along the (0K0), (H20), and (1K0) directions to the numerical evaluation of Eqs. (3) and (5) convolved with a Gaussian approximation to the instrumental energy and wave vector resolution functions while also including the spherical approximation for the Yb³⁺ magnetic form factor and an additive background term. The resulting spectra are shown in Fig. 2d–f and h–i; the continuum response and the sharp feature within the continuum are reproduced exceptionally well. (see  Supplementary material for additional comparisons between the data and the model.) We remark that the excellent agreement between the experimental data and the ideal HLHM is also in accordance with the prediction of ref. ⁶⁴.

The sharp feature toward the top of the continuum is a particularly interesting aspect of the spectrum that, to the best of our knowledge, has not been previously observed in a quantum magnet. In the model, such a sharp feature appears within the two-magnon continuum due to a VHS in the joint density of states^{49,50,52–54}. Indeed, on the level of pure kinematics (i.e., ignoring any matrix element effects), the longitudinal two-magnon response in Eq. (5) is proportional to the joint density of states, $$ {\widehat{g}}_{\mathbf{q}}\left(\omega \right)={\sum }_{\mathbf{k}}\delta \left(\omega -{\epsilon }_{\mathbf{k}}-{\epsilon }_{\mathbf{q}-\mathbf{k}}\right)$$, at each momentum q, which corresponds to the joint band dispersion $$ {\widehat{\epsilon }}_{\mathbf{q}}\left(\mathbf{k}\right)={\epsilon }_{\mathbf{k}}+{\epsilon }_{\mathbf{q}-\mathbf{k}}$$ as a function of the individual magnon momentum k. Being a two-dimensional band dispersion, $$ {\widehat{\epsilon }}_{\mathbf{q}}\left(\mathbf{k}\right)$$ has a saddle-point VHS which gives a logarithmic divergence in the density of states $$ {\widehat{g}}_{\mathbf{q}}\left(\omega \right)$$ and, thus, in the longitudinal spin response. Physically, the VHS is a specific energy transfer ω that can create many distinct magnon pairs with a fixed total momentum q but different individual momenta k and q − k. The coalescence of such distinct scattering processes is analogous to the coalescence of light rays giving rise to caustic features in ray optics (see Fig. 4). Therefore, using this analogy, the VHS can also be understood as a caustic feature in the longitudinal spin dynamics. We emphasize that the observation of a VHS is direct evidence for strong quantum fluctuations in YbCl₃ (because the VHS appears in the longitudinal spin response) as well as the two-dimensional nature of its quantum magnetism (because significant interlayer exchange would smear out the VHS).

Fig. 4

Alternative understanding of Van Hove singularities as caustic features in the collective spin dynamics.

a Caustic features in ray optics. Parallel light rays (black lines) enter an optical system at different positions. When these light rays reflecting from a circular mirror (red line) coalesce, they give rise to caustic features in real space. b Caustic features in a spin response. The two-magnon continuum can be understood as a sum of sharp contributions, ω = ε_k + ε_q−k, each corresponding to a fixed momentum k of the first magnon. When these sharp contributions (black lines) coalesce, they give rise to caustic features in the two-magnon continuum. Note that the spin response shown here is for a one-dimensional model system; for the two-dimensional system in consideration, the caustic features appear inside the continuum (not at its edge) and are weaker as they correspond to logarithmic (rather than square-root) singularities.

Upon close examination, the analytic model does not fully capture the intensity and dispersion of the VHS over the entire zone, as can be seen in Fig. 2a–f. By plotting the fitted peak positions of the VHS from Fig. 2g on the calculated spectra in Fig. 2d–f, we notice that there are differences between the calculated and the observed VHS energies near the (100) and $$ \left(\frac{1}{2}20\right)$$ wave vectors  (≈0.2 meV). These energy differences may arise from a renormalization of the two-magnon continuum due to magnon interactions neglected in linear spin-wave theory. Also, we point out that domain formation due to bond disorder⁶⁵ could, in principle, play a role in partially smearing out the VHS and thus explaining why certain portions of the predicted VHS are absent from the experimental data.

We now explore the potential importance of additional exchange interactions through comparison of the experimental spin-wave dispersion, including points measured along the L direction, to calculations using the SpinW software⁷³. Since there are very small differences in the bond lengths within the honeycomb layers of YbCl₃, as described in Fig. 1, we label two of the three nearest-neighbor exchanges as J for the d = 3.884 Å bonds and the third one as $$ J^{\prime }$$ for the d = 3.867 Å bond. The numerical comparison yields J = 0.42(1) meV and $$ J^{\prime }=0.43\left(1\right)$$ meV, while J_c refines to a value below the detection limit of 0.016 meV for the measurements reported here (see  Supplementary material for additional details). The resulting cross-section is shown in Fig. 3g–k. The numerically determined J and $$ J^{\prime }$$ are indistinguishable from each other and close to the value determined by a comparison to the analytical model. The spin-wave modes from linear spin-wave theory accurately reproduce the dispersion and intensity distribution of this portion of the measured spectrum (Fig. 3g–k). We also attempted to include next-nearest-neighbor exchange interactions within the plane of the honeycomb lattice; the best fit values of J₂ and $$ J_2^{\prime }$$ are three orders of magnitude smaller than J and zero within the error bars of the refined value (see Supplementary material for additional details).

Magnetic moment and Néel temperature

To further check the validity of the ideal HLHM, we have performed polarized neutron diffraction measurements at T =  0.3 K (see Supplementary material for additional details). We observe an antiferromagnetic spin arrangement with a well-defined staggered magnetic moment. By symmetry, this moment can either lie in the ac plane or point along the b axis, i.e., a linear combination of the two directions is not allowed. Our measurements reveal that the magnitude of the ordered moment is 1.06(4) μ_B and that the moment points primarily along the a axis with only a small deviation of 5(3)° along the c axis. The expected moment from the ground state crystal field doublet is 2.24(5) μ_B⁶¹. The value of the ordered moment is thus  ≈47% of the fully polarized moment, which is in excellent agreement with the value of  ≈48% obtained from linear spin-wave theory and the value of  ≈54% determined by more accurate techniques^72,74–77. We emphasize, however, that these theoretical values are extremely sensitive to any interlayer exchange J_c. For example, if we include ∣J_c/J∣ ≈ 4%, corresponding to the detection limit from the INS response (see above), the ordered moment in linear spin wave theory already increases from  ≈48% to  ≈63%. Therefore, the small value of the experimental moment indicates that the actual value of J_c is significantly below the INS detection limit and further confirms our conclusion that the magnetism in YbCl₃ is two dimensional to a very good approximation.

The same conclusion is demonstrated even more strikingly by the very small value of the Néel temperature, T_N ≈  0.6 K⁶². Theoretically, for a quasi-two-dimensional Heisenberg magnet, the Néel temperature is approximately given by $$ T_N=4\pi {\rho }_s/\log \mid J/J_c\mid $$, where ρ_s is the spin stiffness of the purely two-dimensional system. At this temperature, T = T_N, the individual interlayer exchanges J_c add constructively within an area of linear size $$ \xi ~\exp \left(2\pi {\rho }_s/T\right)$$^78,79 to produce an effective interlayer exchange on the order of the intralayer exchange J. Taking ρ_s ≈ 0.101J from the literature^76,77 and J ≈ 0.42 meV from the neutron scattering data, the relative interlayer exchange is then estimated to be $$ \mid J_c/J\mid ~\exp \left(-4\pi {\rho }_s/T_N\right)\approx 3\times 10^{-5}$$. This value is three orders of magnitude smaller than the detection limit from the INS response and indicates that the magnetism in YbCl₃ is as close to two-dimensional as in the canonical cuprate antiferromagnet La₂CuO₄⁸⁰. Moreover, the specific ordering pattern below T_N signals the presence of very small anisotropic interactions on top of the dominant Heisenberg exchange J. Since these anisotropic spin interactions should also contribute to T_N, the interlayer exchange may be even smaller than our estimate. We finally remark that the application of a small magnetic field should, in principle, cause the ordering transition to sharpen and the Néel temperature to increase⁷⁹. While there are experimental indications for such behavior in the specific heat⁶², more detailed studies of the field dependence would be very important to confirm not only this prediction but also further predictions of interesting quantum effects at larger magnetic fields⁵⁴.

Heat capacity calculations

Specific heat capacity measurements provide an additional means of examining the HLHM in YbCl₃. The experimental specific heat and the entropy of YbCl₃ for 0.5 K < T < 8 K are shown in Fig. 5a and b. Between T = 0.5 K and T = 8 K, nearly all of the entropy, $$ R\ln \left(2\right)$$, for the ground state doublet has been recovered by the system with only a very small contribution in the region of the transition to long-range magnetic order⁶². We use mTPQ calculations⁸¹, as implemented in the $$ \mathcal{H}\mathrm{\Phi }$$ library⁸², for a cluster size of 32 spin $$ \frac{1}{2}$$ elements to calculate the heat capacity as a function of the reduced temperature T/J (see Supplementary material for additional details). The results for the HLHM with J = 0.42 meV, obtained by fitting the INS data, are shown in Fig. 5a and b. The overall shape is in reasonable agreement with the data, but a somewhat improved comparison is found by using J = 0.32 meV. The discrepancy may be due to quantum corrections, such as renormalizations due to magnon interactions^51,52, that are captured by the mTPQ calculations but neglected in linear spin wave theory.

Fig. 5

Magnetic specific heat capacity, entropy, and mTPQ calculations for YbCl₃.

Specific heat capacity (a) and entropy (b) for YbCl₃. The solid red line is the best fit calculation using mTPQ as described in the text. The dashed blue line is the calculated specific heat using mTPQ with the value J = 0.42 meV obtained by fitting the INS data.

We have used INS to investigate the collective magnetic excitation spectrum of YbCl₃. In addition to a conventional transverse spin-wave (single-magnon) mode, there is a longitudinal two-magnon continuum harboring a sharp VHS. These components are all reproduced by linear spin-wave theory with a single nearest-neighbor Heisenberg interaction on the honeycomb lattice. The sharp VHS is observable due to the almost perfectly two-dimensional quantum magnetism in YbCl₃, which is further reflected in a strongly reduced ordered moment and a very small ordering temperature. Our results establish YbCl₃ as an ideal model system to investigate the collective quantum behavior of the honeycomb antiferromagnet in the unfrustrated limit. Finally, we point out that the observation of a VHS in a two-magnon continuum provides a strong indication that a similar observation in the two-spinon spectrum of a two-dimensional quantum spin liquid^83,84 is experimentally feasible. Such an observation in a quantum spin liquid would be important in ruling out competing sources of a continuum response, such as quenched disorder or overdamped magnons, and could provide a more unambiguous signature of a long-range-entangled quantum state.