Fractional excitations in the square lattice quantum antiferromagnet.

Quantum magnets have occupied the fertile ground between many-body theory and low-temperature experiments on real materials since the early days of quantum mechanics. However, our understanding of even deceptively simple systems of interacting spins-1/2 is far from complete. The quantum square-lattice Heisenberg antiferromagnet (QSLHAF), for example, exhibits a striking anomaly of hitherto unknown origin in its magnetic excitation spectrum. This quantum effect manifests itself for excitations propagating with the specific wave vector (π, 0). We use polarized neutron spectroscopy to fully characterize the magnetic fluctuations in the metal-organic compound CFTD, a known realization of the QSLHAF model. Our experiments reveal an isotropic excitation continuum at the anomaly, which we analyse theoretically using Gutzwiller-projected trial wavefunctions. The excitation continuum is accounted for by the existence of spatially-extended pairs of fractional S=1/2 quasiparticles, 2D analogues of 1D spinons. Away from the anomalous wave vector, these fractional excitations are bound and form conventional magnons. Our results establish the existence of fractional quasiparticles in the high-energy spectrum of a quasi-two-dimensional antiferromagnet, even in the absence of frustration.

A fascinating manifestation of quantum mechanics is the emergence of elementary excitations carrying fractional quantum numbers. Fractional excitations were a central ingredient to understand the fractional quantum Hall effect [1], and have been investigated in a range of systems including conducting polymers [2], bilayer graphene [3], cold atomic gases [4], and low-dimensional quantum magnets [5, 6]. Among the latter class of systems, the spin-1/2 Heisenberg antiferromagnet chain (HAFC) is perhaps the simplest model, for which the ground-state and the excitations are known exactly [7-9]. Excitations of the spin-1/2 HAFC created by an elementary ΔS = 1 process are radically different from spin-waves, the coherent propagation of a flipped spin, and are pairs of unbound fractional quasiparticles known as spinons, each carrying a S = 1/2 quantum number. The existence of spinons in the spin-1/2 HAFC has been confirmed experimentally in a number of quasi-1D materials [10, 11], but observing their 2D and 3D analogues is an ongoing challenge [6]. To date the main candidate systems comprise geometrically frustrated magnets on the triangular [12] or kagome [13-15] lattices. In this work, we take a frustration-free route and focus on the quantum (spin-1/2) square-lattice Heisenberg antiferromagnet (QSLHAF), one of the most fundamental models in magnetism. It is defined by the Hamiltonian where the antiferromagnetic exchange J is restricted to nearest-neighbor spin pairs 〈i, j〉. We provide experimental and theoretical evidence that even in this simplest of 2D models deconfined fractional S = 1/2 quasiparticles can be identified at high energies, where they modify the short wave-length spin dynamics and are responsible for a significant quantum anomaly that cannot be captured by conventional spin-wave theory.

It may seem surprising that the QSLHAF is a candidate for hosting fractional excitations as at a superficial level its long range magnetic order resembles that of a classical system. The elementary excitations of this ”Néel state”, when calculated using semi-classical spin wave theory (SWT), are bosonic quasiparticles, known as magnons: the one-magnon spectrum is gapless, with two-magnon excitations occupying a continuum at higher energy. The interaction between magnons is relatively weak and leads to an upward renormalization of the magnon energy and to scattering between two-magnon states [16, 17]. One- and two-magnon excitations respectively correspond to fluctuations perpendicular (transverse) and parallel (longitudinal) to the direction of the ordered moments.

While none of the above properties suggest the existence of quasiparticle fractionalization, quantum effects are nevertheless far from negligible in the QSLHAF. This is evidenced by the observation that quantum zero-point fluctuations reduce the staggered moment to only 62% of its fully ordered value S = 1/2 [18, 19]. This suggests that the QSLHAF may in fact be close to a state preserving spin-rotation symmetry, such as the resonating-valence-bond (RVB) state proposed by Anderson [20] for the cuprate realization of this model. In particular, fractional spin excitations present in the RVB state may be relevant for the spin dynamics in the Néel state, especially at high energies. Indeed, analytical theories using bosonic [21] or fermionic [22, 23] fractional quasiparticles have long been proposed and it has been shown that the presence of conventional classical long range order does not hinder the possibility of fractional excitations [24, 25]. By analogy with the 1D case, these are referred to as spinons.

The magnetic excitation spectrum of various realizations of the QSLHAF have been investigated using neutron spectroscopy, including the parent compounds of the high-T cuprate superconductors Sr2CuO2Cl2 [26, 27] and La2CuO4 [28, 29], Sr2Cu3O4Cl2 [30] and the metal-organic compounds Cu(pz)2(ClO4)2 [31, 32] and Cu(DCOO)2·4D2O (CFTD) [33, 34] considered here. These experiments have established that while SWT gives an excellent account of the low-energy spectrum, a glaring anomaly is present at high energy for wavevectors in the vicinity of (π, 0), where q = (q) is expressed in the square-lattice Brillouin-zone of unit-length 2π. The anomaly is evident as a dramatic wipe out of intensity (Fig. 1a) of the otherwise sharp excitations [27, 29, 32, 34] and as a 7% downward dispersion along the magnetic zone-boundary connecting the (π/2, π/2) and (π, 0) wavevectors for Sr2Cu3O4Cl2 [30, 33] and CFTD (see also Fig. 4d). Unambiguously identifying the origin of this effect is complicated by the presence, in some of these materials, of small additional exchange terms such as electronic ring-exchange [27, 29], further neighbour exchange [31, 32] or interpenetrating sublattices [30]. In contrast, the deviations observed in CFTD agree with numerical results obtained by series expansion [35, 36], quantum Monte Carlo [37, 38], and exact diagonalization [39] methods for the model of Eq. (1), proving that the anomaly is in this case intrinsic [34]. Due to the similarities of the measured anomaly with some aspects of the predicted fermionic RVB excitations treated in the random phase approximation [23], it has been speculated that the anomaly might be related to fractionalized spin excitations [29, 34]. Given the dramatically enlarged family of experimentally accessible QSLHAF physical realization due to the advent of high-resolution Resonant Inelastic X-ray Scattering [40-43] and the fundamental nature of the QSLHAF, it is clearly desirable to develop a microscopic understanding of the origin of the anomaly.

Figure 1

Overview of the magnetic excitation spectrum of CFTD and its interpretation in terms of spin-waves or spatially-extended fractional excitations

(a) Momentum and energy dependence of the (total) dynamic structure factor measured by time-of-flight inelastic neutron scattering. Square boxes (black dashed) highlight the (π, 0) and (π/2, π/2) wave vectors. (b) and (d) Corresponding distributions of real-space fractional quasiparticle-pair separations, as calculated in the ∣SF〉 variational state (Eq. 6), evidencing respectively the unbound and bound nature of the pair excitations. (c and e) Pictorial representation of a quasiparticle pair excitation and a spin-wave excitation (magnon) respectively.

Figure 4

Variational excitation spectra of the Gutzwiller-projected trial wavefunctions

Transverse dynamic structure factor for the ∣SF+N〉 (a) and ∣SF〉 (b) states with lattice sizes of 24 × 24 and 32 × 32 respectively. (c) The magnon-like dispersion extracted from a (red points) compared to the experimental CFTD data [34] (blue squares, error bars corresponds to one standard deviation), spin-wave theory with first-order (solid black line) and third-order [62] (green triangles) 1/S corrections, series expansion [36] (dashed purple line) and quantum Monte Carlo [38] (cyan diamonds). The experimental data is scaled using J = 6.11 meV. (d) Zoom-in on the magnon-like mode dispersion along the magnetic zone-boundary.

Here we present polarized neutron scattering results on CFTD which establish the existence of a spin-isotropic continuum at (π, 0), which contrasts sharply with the dominantly longitudinal continuum at (π/2, π/2) and with the broken spin-symmetry of the ground state. Using a fermionic description of the spin dynamics based on a Gutzwiller-projected variational approach, we argue that the continuum at (π, 0) is a signature of spatially extended pairs of fractional S = 1/2 quasiparticles (Fig. 1b and 1c). At other wave vectors, including (π/2, π/2) (Fig. 1d), our approach yields bound pairs of these fractional quasiparticles and so recovers a conventional magnon spectrum, in agreement with SWT (Fig. 1e).

Neutron scattering experiments were performed on single crystals of CFTD using unpolarized time-of-flight spectroscopy (Fig. 1) and triple-axis spectroscopy with longitudinal polarization analysis (see Supplementary Materials). The results of our polarized experiment are presented in Fig. 2 through the energy dependence of the diagonal components of the dynamic structure factor. By combining wave vectors from different equivalent Brillouin zones (see Supplementary Materials), we can reconstruct the total dynamic structure factor (Fig. 2a and e) and separate contributions from spin fluctuations that are transverse to and along (Fig. 2b-c and f-e) the ordered moment. Within SWT, the resulting transverse and longitudinal spectra are dominated by one-magnon and two-magnon excitations, respectively. At (π/2, π/2) and at an energy of ω/J = 2.38(2) we observe a sharp, energy resolution limited peak (Δω=1.47(5) meV=0.24(1) J, FWHM) which is the signature of a long-lived, single-particle excitation (Fig. 2e). Most of the observed spectral weight is in the resolution-limited peak of the transverse channel (Fig. 2f), while a weak continuum extends from ω/J ≈ 2.3 to 3.4 with a maximum around ω/J ≈ 2.6 in the longitudinal channel, (Fig. 2g). In contrast, the response at (π, 0) displays a pronounced high-energy tail, starting right above the peak maximum at ω/J = 2.19(2) and extending up to ω/J ≈ 3.8. This tail carries 40(12)% of the total spectral weight at (π, 0) (Fig. 2a) and is evident in both the transverse (Fig. 2b) and longitudinal (Fig. 2c) channels. To isolate the continuous component in the transverse channel we subtract resolution-limited Gaussians corresponding to sharp, single-particle responses, with the results shown in Fig. 2d and h. This analysis reveals the important fact that the transverse continuum at (π, 0) is within error twice the longitudinal contribution (Fig. 2d). Thus we can conclude that the continuum at (π, 0) arises from correlations which are isotropic in spin space with , while by contrast the continuum contribution at (π/2, π/2) is fully contained in the longitudinal channel (Fig. 2h).

Figure 2

Summary of the polarized neutron scattering data

Energy dependence of the total, transverse and longitudinal contributions to the dynamic structure factor at constant wave vectors (a–c) q = (π, 0) and (e–g) q = (π/2, π/2) measured by polarized neutron scattering on CFTD. The solid lines indicate resolution-limited Gaussian fits, while the dashed lines are empirical lineshapes used as guides-to-the-eye. (d and h) transverse dynamic structure factor with subtracted resolution-limited Gaussian fits at (π, 0) and (π/2, π/2) respectively. Error bars correspond to one standard deviation

The pronounced asymmetric and non-Lorentzian line shape of the continuum at (π, 0) cannot be accounted for by conventional effects, even including instrumental resolution. SWT predicts that magnon interactions transfer up to 20% of the transverse spectral weight at the zone-boundary from the sharp one-magnon peak to a higher energy continuum of three-magnon states [17]. However, the resulting line shape differs radically from our observations, does not coincide with the longitudinal response, and does not appear to depend significantly on the wave vector along the zone boundary. Spontaneous magnon decays can in principle produce an asymmetric line shape but are prohibited in this case by the collinearity of the magnetic order [16, 44]. Instead, recent quantum Monte Carlo work [45] suggests to look for explanations of the continuum contribution to the dynamic structure factor at (π, 0) involving the deconfinement of fractional excitations. This is further motivated by the observed coexistence of sharp two-spinon bound states with a broad multi-spinon continuum, at comparable energy ranges but different wave vectors, in the quasi-2D materials Cs2CuCl4 [12, 46] and LiCuVO4 [47] made of strongly-coupled Heisenberg chains.

In order to explore whether fractionalization of magnons can account for the (π, 0) anomaly in the QSLHAF, we use a theoretical approach based on Gutzwiller-projected variational wave functions [48, 49]. In this approach, spin operators are transformed into pairs of S = 1/2 fermionic operators so that Eq. 1 becomes where creates (annihilates) an electron with spin σ at site i. This transformation embeds the original spin Hilbert space into an electronic Hilbert space which also contains nonmagnetic sites occupied by zero or two electrons. As a result, Eqs 1 and 2 are only equivalent on the restricted electronic subspace with half electron filling and no empty sites or double occupancies. This constraint can be enforced exactly by the so-called Gutzwiller projector P. The advantage of this approach is that pairs of fractional S = 1/2 quasiparticles (for the original spin Hamiltonian) can be naturally constructed as particle-hole excitations in the electronic space, projected a-posteriori by P onto spin configurations with exactly one electron per site [50]. The projection may be approximated using the Gutzwiller approximation [22] or the random phase approximation [23]. In this work we choose to implement the projection exactly using the numerical Variational Monte Carlo technique [49, 51].

The quartic electronic operator in Eq. 2 is treated by a mean-field decoupling where the averages and are considered. We adopt the following Ansätze for their real-space dependencies : corresponds to a staggered Néel order parameter (N) and to a staggered flux (SF) threading square plaquettes of the lattice [52-54] (see Supplementary Materials for exact definitions and more details). To each average corresponds a variational parameter whose value is optimized to minimize the energy (Eq. 1) of the Gutzwiller-projected state, ∣SF+N〉 = P∣ψMF〉. The corresponding mean-field electronic ground-state ∣ψMF〉 contains empty and doubly occupied sites and reads where ∣0〉 is the electron vacuum and where the operators are linear combinations of operators that diagonalize the mean-field electronic Hamiltonian. The product over the wave vector k is restricted to the magnetic Brillouin zone (MBZ), a result of the antiferromagnetic unit-cell-doubling. Consequently “±” denotes the band index. In the present case of half electron filling, the “−” band is fully occupied, and there is a finite gap to the empty “+” band for non-zero Néel order-parameter. The overall minimization procedure is carried out numerically using Variational Monte Carlo and leads to a ∣SF+N〉 state with variational energy ESF+N = −0.664J and staggered moment 0.75S per site [48, 55]. This can be compared to more precise Green’s function Monte Carlo studies for Eq. 1 that obtained −0.669J and 0.615S for the ground-state energy and the staggered moment, respectively [56, 57].

The optimized ∣SF+N〉 state, while giving a good estimate for the ground-state energy, does not have the correct long-distance behaviour for the transverse equal-time correlator 〈S+(0)S−(r)〉, predicted by SWT to decay as a power-law [16]. This algebraic decay is a robust long-wavelength prediction and has been implemented in variational magnetic trial wavefunctions in the past [58, 59]. Instead, as the excitation spectrum of the mean-field electronic ground-state is gapped, 〈S+(0)S−(r)〉 decays exponentially after projection. We conjecture that the asymptotic behaviour of the spin correlator is important for the deconfinement of fractional excitations. To obtain insight into the influence of long-distance spin fluctuations, we consider a distinct variational state, ∣SF〉, for which the finite staggered-flux is retained but the Néel order is reduced to zero. ∣SF〉 is a quantum spin-liquid singlet of variational energy ESF = −0.638J that displays a power-law decay of its transverse spin correlations [60, 61].

We now turn to the construction of transverse (S = 1) spin excitations for the above variational states, aiming at comparing their respective dynamic structure factor with the results of Fig. 2. The variational transverse spin excitations are obtained as superpositions of projected particle-hole excitations with momentum q: where the states ∣k, q〉 are generated by destroying a spin-down quasiparticle in the “−” band and creating a spin-up quasiparticle in the “+” band. The coefficients are obtained by diagonalizing the original Hamiltonian (Eq. 1) projected onto the non-orthonormal set of states ∣k, q〉 and correspond to the eigen-energies . Expressing the Fourier-space quasiparticle operators γk using the real-space c operators, we note that the variational spin excitations contain both local spin flips (Fig. 3b) and spatially separated particle-hole excitations, (Fig. 3c). The dynamic structure factor of the transverse spin excitations is calculated as where ∣GS〉 stands either for ∣SF+N〉 or ∣SF〉. We use the identity valid for both variational ground-states to compare the transverse dynamic structure factor of the variational states ∣SF+N〉 and ∣SF〉 with the experimental results presented in Fig. 2. A similar approach also allows to obtain the longitudinal (S = 0) dynamic structure factor (see Supplementary Materials).

Figure 3

A schematic representation of local-spin-flip and spatially separated quasiparticle pair excitations in the Gutzwiller-projected approach

(a) The mean-field wave function ∣ψMF〉 is shown as a resonating-valence-bond liquid (for a better visualization, all singlets are shown as nearest-neighbour and the Néel order is ignored). Configurations containing doubly occupied sites (right hand side) are discarded by the Gutzwiller projection P. (b) Local spin flips create triplets out of resonating singlets. Configurations from ∣ψMF〉 originally containing doubly occupied sites are still projected out (right hand side). (c) Non-local quasiparticle-pair excitations are constructed as projected particle-hole excitations. At a non-zero separation r, they contribute by annihilating a doubly occupied site with a hole, leaving two separated spins-↑. After projection, the only configurations left are the ones constructed from ∣ψMF that contained one empty and one doubly occupied site (right hand side).

The transverse dynamic structure factor of the ∣SF+N〉 state is shown in Fig. 4(a), as obtained by variational Monte Carlo on a finite lattice of 24 × 24 sites. The dominant feature of the spectrum is a low-energy magnon-like mode, which resembles the experimental results of Fig. 1a. In particular, our calculation produces a dispersion along the magnetic zone-boundary in better quantitative agreement with the 7% dispersion observed in Ref. [34] than any other theoretical method, Figs. 4c and 4d. This confirms that magnons can be quantitatively interpreted as bound pairs of fractional S = 1/2 quasiparticles.

However the ∣SF+N〉 transverse dynamic structure factor exhibits a gap at (π, π) and no continuum above the magnon branch at (π, 0). We believe that this is an artefact of replacing the spontaneous symmetry breaking by a Néel mean-field order parameter: this ansatz, as mentioned above, distorts the long-distance asymptotics of spin correlations. Indeed, if we reduce the Néel mean-field parameter of the ∣SF+N〉 state, then the high-energy excitations at (π, 0) move down in energy (see Supplementary Materials). When the Néel field vanishes, (i.e., in the ∣SF〉 state) they evolve into a succession of modes distributed on an extended energy range above the spin-wave mode (shown in Fig. 4b for a 32 × 32 lattice). This behaviour contrasts the situation at (π/2, π/2) where the high-energy transverse excitations completely lose their spectral weight on reducing the Néel field and only the spin-wave mode remains in the ∣SF〉 state. At (π, π), the lowest mode moves down reaching negative energy, which indicates an instability of the ∣SF〉 state towards Néel ordering. We therefore suggest that the continuum of excitations observed at (π, 0) is conditionally dependent on power-law transverse spin correlations and that it corresponds to deconfined fractional spin-1/2 quasiparticles.

To support this interpretation, we consider in Fig. 5a and 5b the system-size dependence of . While the excitations at (π/2, π/2) form a single sharp mode with energy and intensity nearly independent of the system size, the number of modes at (π, 0) and their relative weights are strongly modified by increasing the number of sites. This behaviour is consistent with the development of a continuum of fractional quasiparticles at (π, 0) in the thermodynamic limit.

Figure 5

Finite-size effects and real-space structure in the ∣SF〉 state

Transverse dynamic spin structure factor at (π, 0) (a) and (π/2, π/2) (b) for different system sizes ranging from 8 × 8 (dark blue line) to 32 × 32 (dark red line). Disk-integrated fractional-quasiparticle pair separation distribution Pq(r) at (π, 0) (c) and (π/2, π/2) (d) for corresponding system sizes. (e) Mean fractional-quasiparticle pair separation rq at (π, 0) (red symbols) and (π/2, π/2) (blue symbols). Error bars correspond to one standard deviation from the variational Monte Carlo sampling.

Having established that our Gutzwiller approach depending on wave vector produces respectively sharp and continuum-like excitations from the ∣SF〉 state, we analyze their real-space structure to gain further insight into their nature. We consider their overlap with projected real-space particle-hole excitations , where a Fourier transformation was applied to reflect translation invariance. In this formalism, the most local projected particle-hole pair is the spin-flip state corresponding to a magnon while non-local pairs are characterized by a finite separation r. Therefore, the degree of deconfinement of a fractional S = 1/2 quasiparticles pair can be characterized using the spatial extent of its overlap with projected real-space particle-hole excitations 〈q, r, +∣q, n, +〉. Since the continuum in 5a is populated by different sets of discrete modes for the various system sizes considered, we choose to evaluate the degree of deconfinement through a single q-specific averaged quantity where the aforementioned overlap is weighted by the intensity of the transverse spin excitation in the dynamic structure factor, thus only accounting for modes proportionally to their spectral weight. The spatial profile of ρq(r) for the magnetic zone-boundary wave vectors, shown in Figs. 1b and 1d, reveals much more extended fractional S = 1/2 quasiparticles pairs at (π, 0) than at (π/2, π/2). This is confirmed by the system-size dependence of the radially-integrated normalized distribution Pq(r) = ∑∣r’∣<ρq(r’), plotted in Figs. 5c and 5d. At (π/2, π/2), Pq(r) saturates at a distance, r, that is nearly independent of the system size, while at (π, 0) it does so at a distance that increases with the number or sites. Similarly, the “root-mean-square” fractional quasiparticles pair separation rq = [∑r∣r∣2ρq(r)/∑rρq(r)]1/2, presented in Fig. 5e, grows nearly linearly with the system size for (π, 0) while it has a much weaker size dependence for (π/2, π/2).

Taken together, our real-space results for the ∣SF〉 state show that spin excitations at (π/2, π/2) can indeed be considered as bound pairs of S = 1/2 quasiparticles with confined spatial extent. In contrast at (π, 0) the strong system size dependence of the spin excitations spatial extent indicates the perhaps only marginal deconfinement of fractional quasiparticles in two spatial dimensions. Note that even in the absence of long range Néel order, deconfinement only happen at the special point (π, 0) and no continuum develops at (π, π) as would be naively expected for an algebraic spin liquid. This indicates that the deconfined (π, 0) excitations should only be considered as remnants of the underlying unprojected deconfined particle-hole excitations. This suggests that the QSLHAF ground state can still be understood as a conventional Néel state different from the AF* state described in Refs. 24, 25, where magnons and spinons represent two different branches of excitations. We do not attempt to extract power laws from the numerical data, since the variational ∣SF〉 state mimics the long-distance spin correlations only qualitatively.

Combining our polarized neutron scattering and theoretical results provides evidence that even in the simplest of 2D spin models, deconfined fractional S = 1/2 quasiparticles can be identified at high energies, and account for the quantum anomaly observed in a broad range of experimental realizations of the square-lattice Heisenberg antiferromagnet. This insight raises important theoretical and experimental questions. First, how to obtain explicit quasiparticle deconfinement out of the magnetically ordered ground state of the QSLHAF? How will the excitations uncovered here evolve upon weakening magnetic exchange in one direction hence approaching the 1D limit? Our present work focused on the nearest-neighbor Heisenberg model, an insulator obtained in the strong Coulomb repulsion limit of a one-band Hubbard model. It will be interesting to track the fractional quasiparticles in systems closer to an insulator-metal transition and eventually upon doping. Given that fractional spin excitations are identified at high energies, one may speculate whether weak 2D Mott insulators could, in certain areas of momentum space, host a phenomenon similar to the observed spin-charge separation in 1D [63].