Nearly Deconfined Spinon Excitations in the Square-Lattice Spin-$1/2$ Heisenberg Antiferromagnet

We study the spin-excitation spectrum (dynamic structure factor) of the spin-$1/2$ square-lattice Heisenberg antiferromagnet and an extended model (the $J\text{−}Q$ model) including four-spin interactions $Q$ in addition to the Heisenberg exchange $J$. Using an improved method for stochastic analytic continuation of imaginary-time correlation functions computed with quantum Monte Carlo simulations, we can treat the sharp ($\delta$-function) contribution to the structure factor expected from spin-wave (magnon) excitations, in addition to resolving a continuum above the magnon energy. Spectra for the Heisenberg model are in excellent agreement with recent neutron-scattering experiments on $\mathrm{Cu}(\mathrm{DCOO}{)}_2·4{\mathrm{D}}_2\mathrm{O}$, where a broad spectral-weight continuum at wave vector $q=(\pi ,0)$ was interpreted as deconfined spinons, i.e., fractional excitations carrying half of the spin of a magnon. Our results at $(\pi ,0)$ show a similar reduction of the magnon weight and a large continuum, while the continuum is much smaller at $q=(\pi /2,\pi /2)$ (as also seen experimentally). We further investigate the reasons for the small magnon weight at $(\pi ,0)$ and the nature of the corresponding excitation by studying the evolution of the spectral functions in the $J\text{−}Q$ model. Upon turning on the $Q$ interaction, we observe a rapid reduction of the magnon weight to zero, well before the system undergoes a deconfined quantum phase transition into a nonmagnetic spontaneously dimerized state. Based on these results, we reinterpret the picture of deconfined spinons at $(\pi ,0)$ in the experiments as nearly deconfined spinons—a precursor to deconfined quantum criticality. To further elucidate the picture of a fragile $(\pi ,0)$-magnon pole in the Heisenberg model and its depletion in the $J\text{−}Q$ model, we introduce an effective model of the excitations in which a magnon can split into two spinons that do not separate but fluctuate in and out of the magnon space (in analogy to the resonance between a photon and a particle-hole pair in the exciton-polariton problem). The model can reproduce the reduction of magnon weight and lowered excitation energy at $(\pi ,0)$ in the Heisenberg model, as well as the energy maximum and smaller continuum at $(\pi /2,\pi /2)$. It can also account for the rapid loss of the $(\pi ,0)$ magnon with increasing $Q$ and the remarkable persistence of a large magnon pole at $q=(\pi /2,\pi /2)$ even at the deconfined critical point. The fragility of the magnons close to $(\pi ,0)$ in the Heisenberg model suggests that various interactions that likely are important in many materials—e.g., longer-range pair exchange, ring exchange, and spin-phonon interactions—may also destroy these magnons and lead to even stronger spinon signatures than in $\mathrm{Cu}(\mathrm{DCOO}{)}_2·4{\mathrm{D}}_2\mathrm{O}$.

Popular Summary

Similar to how sound waves propagate through everyday materials, some materials can carry magnetic waves, which are collective excitations of electronic magnetic moments (or spins). These “spin waves” can be detected in neutron-scattering experiments. Theoretically, it is believed that the spin waves can sometimes split up (or fractionalize) into a different kind of excitation, called a spinon, which carries only half of the magnetic moment associated with a normal spin wave. Recent experiments have revealed anomalous spin-wave behavior in antiferromagnetic materials, where the electron spins are oriented antiparallel to their neighbors in the crystal lattice, like the colors of a checkerboard. This behavior suggests that spin waves propagating in certain directions are unstable and fractionalize into spinons. We demonstrate that this proposal is only partially correct.

Using a theoretical quantum spin model known to be an accurate representation of the experimental system and a novel computer-simulation method for extracting the dynamical signatures, we show that the anomalous spin waves are not quite broken up into spinons but exhibit a “mixed personality,” fluctuating between spinons and spin waves. However, by changing the parameters of the model, we can cause the spin wave to fully fractionalize. We also present a simplified theoretical picture in which the mechanism of fractionalization can be captured and understood.

Our achievement demonstrates that spinons are more common than previously thought. We argue that further signs of spin-wave fragility are present in existing neutron-scattering data, and we suggest experiments to test our proposed fractionalization mechanism.

Figure 1

Parametrizations of the spectral function used in this work. In diagram (a), a large number of $\delta$ functions with the same amplitude occupy frequencies ${\omega }_i$ in the continuum (or, in practice, on a very fine frequency grid). The locations are sampled in the SAC procedure. In diagram (b), the $\delta$ function at the lowest frequency ${\omega }_0$ has a larger amplitude, $a_0>a_i$ for $i>0$, and this amplitude is optimized in the way described in the text. The frequencies of all the $\delta$ functions, including ${\omega }_0$, are sampled as in (a) but with the constraint ${\omega }_0<{\omega }_i\forall i>0$.

Figure 2

The goodness of fit vs the amplitude $a_0$ of the lowest $\delta$ function in three runs with different noise realizations for a synthetic spectrum with a $\delta$ function of weight $a_0=0.4$ at ${\omega }_0=1$. The continuum is a Gaussian of width 1 centered at the same ${\omega }_0$, with the weight below ${\omega }_0$ excluded. The noise level is ${\sigma }_i\approx {10}^{-5}$, and the errors are correlated with autocorrelation time 1 according to the description in the Appendix. The inset shows the data close to the $⟨{\chi }^2⟩$ minimum on a different scale.

Figure 3

Mean value of ${\omega }_0$ in SAC runs with four different noise realizations (shown in different colors) graphed vs the amplitude parameter $a_0$. The noise levels are ${10}^{-5}$ and ${10}^{-6}$ in panels (a) and (b), respectively, and the number of $\delta$ functions was $N_{\omega }=1000$ and $N_{\omega }=2000$. In both cases, $G(\tau )$ values were generated on a uniform grid with ${\mathrm{\Delta }}_{\tau }=0.1$ for $\tau$ up to the point where the relative error exceeds 10%.

Figure 4

Two typical SAC-computed spectral functions (red and blue curves, obtained with different noise realizations) compared with the underlying true synthetic spectrum (thicker black curve, with the half-Gaussian containing 60% of the weight). The parameters of the spectrum are the same as in Fig. 6. The noise levels are ${10}^{-5}$ and ${10}^{-6}$ in panels (a) and (b), respectively.

Figure 5

The dynamic structure factor of the 2D Heisenberg model computed on an $L=48$ lattice along a standard path in the BZ indicated on the $x$ axis. The $y$ axis is the energy transfer $\omega$ in units of the coupling $J$. The magnon peak ($\delta$ function) at the lower edge of the spectrum is marked in white irrespective of its weight, while the continuum is shown with color coding on an arbitrary scale where the highest value is 1. The upper white curve corresponds to the location where, for given $\mathbf{q}$, 5% of the spectral weight remains above it.

Figure 6

Dynamic structure factor for the $L=48$ system at four different momenta. The smallest momentum increment $2\pi /L$ is denoted by $k$ in diagrams (a) and (d). The relative amplitude of the magnon pole is indicated in each panel.

Figure 7

Size dependence of the single-magnon energy (a) and weight in the magnon pole (b) at wave vectors $\mathbf{q}=(\pi ,0)$, $(\pi /2,\pi /2)$, and $(\pi ,\pi )$. Lanczos ED results for small systems ($L\times L$ lattices with $L=4$ and $L=6$ as well as tilted lattices with $N=20$, 32, and 40 sites) are shown as open circles, and QMC-SAC data are presented as solid circles with error bars. The error bars were estimated by bootstrap analysis (i.e., carrying out the SAC procedure multiple times with random samples of the QMC data bins).

Figure 8

Comparison of the CFTD experimental data [33] (the full scattering cross section corresponding to unpolarized neutrons) and our QMC-SAC spectral functions at wave vectors $\mathbf{q}=(\pi ,0)$ and $\mathbf{q}=(\pi /2,\pi /2)$. To account for experimental resolution, we have convoluted the QMC-SAC spectral functions in Figs. 6 and 6 with a common Gaussian broadening (half-width $\sigma =0.12J$). We have renormalized the exchange constant by a factor 1.02 relative to the original value in Ref. [33], and to match the arbitrary factor in the experimental data, we have further multiplied both of our spectra by a factor of approximately 50.

Figure 9

Single-magnon dispersion ${\omega }_{\mathbf{q}}$ along a representative path of the BZ. The CFTD experimental data from Ref. [33] are shown as blue squares, and the QMC-SAC data (the location of the magnon pole) are shown with red circles. We also show the linear SWT dispersion (black curve) adjusted by a common factor corresponding to the exact spin-wave velocity $c=1.65847$ [65].

Figure 10

Relative spectral weight of the single-magnon pole along the representative path in the BZ for the $L=48$ Heisenberg system. Error bars were estimated by bootstrapping.

Figure 11

Spectral functions at $\mathbf{q}=(\pi ,0)$ and $(\pi /2,\pi /2)$ obtained using unconstrained SAC with the parametrization in Fig. 1. The insets show comparisons with the experimental data [33], where we have only adjusted a common amplitude to match the areas under the peaks.

Figure 12

Spectral functions obtained using sampling with the parametrization in Fig. 1 under the constraint that no weight falls below the lower bounds determined with a $\delta$ function at the lower edge (Fig. 6): ${\omega }_q=2.13$ and 2.40 for $\mathbf{q}=(\pi ,0)$ and $(\pi /2,\pi /2)$, respectively. The inset shows the results on a different scale to make the continua more visible. The insets of the inset show comparisons with the experimental data, where we have broadened the numerical results by Gaussian convolution and adjusted a common amplitude.

Figure 13

Results for the $J\text{−}Q$ model at $\mathbf{q}=(\pi ,0)$ and $(\pi /2,\pi /2)$, calculated on the $L=32$ lattice. The lowest excitation energy ${\omega }_{\mathbf{q}}$ (a) and the relative weight of the single-magnon contribution (b) are shown as functions of the coupling ratio $Q/J$ from the Heisenberg limit ($Q/J=0$) to the deconfined quantum-critical point ($Q_c/J\approx 22$).

Figure 14

Size dependence of the excitation energy ${\omega }_{\mathbf{q}}$ (a) and the relative weight of the magnon pole $a_0(\mathbf{q})$ (b) at $\mathbf{q}=(\pi ,0)$ close to the Heisenberg limit of the $J\text{−}Q$ model.

Figure 15

The $\mathbf{q}=(\pi ,0)$ dynamic structure factor of the $J\text{−}Q$ model at $Q/J=4$ obtained using SAC with the two parametrizations of the spectrum in Figs. 1 and 1. The relative weight of the leading $\delta$ function in (b) is 1.4%.

Figure 16

(a) Dispersions of the bare excitations of the effective model along a path through the BZ. The lower branch is for the magnon, and the upper branch is for a single spinon. The latter is also the lower edge of the two-spinon continuum. In this example, the spinons in the circled region close to $\mathbf{q}=(\pi ,0)$ almost touch the magnon band, leading to significant spinon-magnon mixing. (b) The black curve shows the lowest energy of the mixed spinon-magnon system obtained with the dispersions in panel (a) and strength $g=5.1$ of the mixing term. The red circles show the results of the QMC-SAC calculations for the Heisenberg model on the $L=48$ lattice from Sec. 3.

Figure 17

Illustration of the mixing process between the magnon (black circle) and the spinon pair (red circles). With mixing strength $g$, a magnon on a given sublattice splits up into a spinon pair occupying nearest-neighbor sites. The spinon pair can recombine and form a magnon on the original sublattice.

Figure 18

Energy levels versus the total wave vector of two spinons interacting through a potential $V(r)=-6.2{\mathrm{e}}^{-r/2}$. The bare dispersion relation of the single spinon is given by Eq. (21) with $c^s=3.1$. We only show a few of the levels between the lower and higher energy bound.

Figure 19

Dispersion relation (a) and wave-vector dependence of the relative weight of the magnon pole (b) calculated with the effective Hamiltonian with the parameters $c^m=3.1$, $c^s=3.1$, $\mathrm{\Delta }=1.94$, $g=1.86$, and the spinon-spinon potential $V(r)=-6.2{\mathrm{e}}^{-r/2}$.

Figure 20

Spectral functions of the effective model at (a) $\mathbf{q}=(\pi ,0)$ and (b) $\mathbf{q}=(\pi /2,\pi /2)$, using model parameters corresponding to the Heisenberg model: $c^m=3.1$, $c^s=3.1$, $\mathrm{\Delta }=1.94$, $g=1.86$, and the spinon-spinon potential $V(r)=-6.2{\mathrm{e}}^{-r/2}$ (same as used in Figs. 18 and 19). The $\delta$ functions in the exact spectral function (computed here using an $L=64$ lattice) have been broadened for visualization.

Figure 21

Spectral functions as in Fig. 20 but with the parameters of the effective model chosen to give behaviors similar to the $J\text{−}Q$ model with $Q\approx J$: $c^m=3.1$, $c^s=6.2$, $\mathrm{\Delta }=0.39$, $g=1.86$, and the spinon-spinon potential $V(r)=-6.2{\mathrm{e}}^{-r/2}$.

Figure 22

(a) Imaginary-time correlation function at $q=(\pi ,0)$ for a 2D Heisenberg lattice with $L=48$, computed in SSE QMC simulations at $\beta =192$ (giving $T=0$ results for all practical purposes). The two straight lines correspond to the contribution from the leading $\delta$ function obtained in the SAC procedure, with amplitude $a_0=0.405\pm 0.025$. (b) The statistical errors (diagonal elements of the covariance matrix) $\sigma (\tau )$ and the eigenvalues of the covariance matrix (ordered from smallest to largest).

Figure 23

The eigenvectors corresponding to the smallest ($n=1$) and largest ($n=40$) eigenvalues of the covariance matrix in Fig. 22, as well as one from the middle of the eigenvalue spectrum ($n=20$).

Figure 24

The same kind of data as in Fig. 22 but obtained using a synthetic spectrum with a $\delta$ function of weight $a_0=0.4$ at ${\omega }_0=1$ and a continuum consisting of a half-Gaussian above the $\delta$ function, of width 1. The straight line in (a) corresponds to the contribution from just the $\delta$ function.