Spinon confinement in a quasi one dimensional anisotropic Heisenberg magnet

TL;DR

This paper investigates spinon confinement in a quasi-one-dimensional antiferromagnetic material using neutron scattering, revealing bound states and validating a new theoretical approach that matches experimental data.

Contribution

It provides the first detailed experimental and theoretical analysis of spinon confinement in SrCo2V2O8, introducing a novel tangent-space matrix product states technique.

Findings

Observation of spinon excitations above TN

Detection of meson-like bound states below TN

Excellent agreement between theory and neutron scattering data

Abstract

Confinement is a process by which particles with fractional quantum numbers bind together to form quasiparticles with integer quantum numbers. The constituent particles are confined by an attractive interaction whose strength increases with increasing particle separation and as a consequence, individual particles are not found in isolation. This phenomenon is well known in particle physics where quarks are confined in baryons and mesons. An analogous phenomenon occurs in certain magnetic insulators; weakly coupled chains of spins S=1/2. The collective excitations in these systems is spinons (S=1/2). At low temperatures weak coupling between chains can induce an attractive interaction between pairs of spinons that increases with their separation and thus leads to confinement. In this paper, we employ inelastic neutron scattering to investigate the spinon confinement in the quasi-1D S=1/2 XXZ antiferromagnet SrCo2V2O8. Spinon excitations are observed above TN in quantitative agreement with established theory. Below TN the pairs of spinons are confined and two sequences of meson-like bound states with longitudinal and transverse polarizations are observed. Several theoretical approaches are used to explain the data. A new theoretical technique based on Tangent-space Matrix Product States gives a very complete description of the data and provides good agreement not only with the energies of the bound modes but also with their intensities. We also successfully explained the effect of temperature on the excitations including the experimentally observed thermally induced resonance between longitudinal modes below TN ,and the transitions between thermally excited spinon states above TN. In summary, our work establishes SrCo2V2O8 as a beautiful paradigm for spinon confinement in a quasi-1D quantum magnet and provides a comprehensive picture of this process.

Equations

$$H$$

$${\omega }_l(Q)=\{\begin{array}{ll}{\omega }_{-}(Q)& \text{for }0\le Q<Q_{\kappa }\\ \frac{2I}{1+\kappa }\sin (Q)& \text{for }{Q}_{\kappa }<Q<\pi /2\\ {\omega }_{1sp}(Q)& \text{for }\pi /2<Q<\pi \end{array}.$$

$${\omega }_u(Q)=\{\begin{array}{ll}{\omega }_{1sp}(Q)& 0<Q<Q_{ϵ}\\ {\omega }_{-}(Q)& Q_{ϵ}<Q<\pi \end{array}$$

$${\omega }_{\pm }(Q)$$

$${\omega }_{1sp}(Q)$$

$$\frac{K(k^{\prime})}{K\big{(}\sqrt{1-k^{\prime 2}}\big{)}}=\frac{1}{\pi}{\rm arccosh}\big{(}\epsilon^{-1}\big{)}.$$

$${\omega }_{\mathrm{V}}(Q)\simeq {\max }_p|{\omega }_{1sp}(p)-{\omega }_{1sp}(p+Q)|.$$

$${\omega }_{\mathrm{V}}(Q)\approx 0.94J\sin (Q).$$

$$I({{\bf Q}^{\prime}},\omega)=\sum_{l=1}^{4}A_{l}({\bf Q}^{\prime})I_{\rm 1D}\left(4Q+\frac{2\pi(l-1)}{c},\omega\right)\$$

$$-\frac{{\hslash }^2}{\mu }\frac{d^2\phi }{d{x}^2}+\lambda |x|\phi =(E-2E_0)\phi .$$

$$E_j=2E_0+\alpha {\zeta }_{j}j=1,2,3,....,$$

$${\Omega }_j^{(a)}(Q)\approx E_j^{(a)}+\frac{(Q-\pi {)}^2}{2m_j^{(a)}},a=L,T,$$

$${\Omega }_{j,k}^{(a\to b)}(Q-P)=|{\Omega }_j^{(a)}(Q)-{\Omega }_j^{(b)}(P)|,a,b=T,L.$$

$$\sum _{i,j,n,m}{J}_{n,m}^{i,j}[S_{i,j}^z{S}_{n,m}^z+ϵ(S_{i,j}^x{S}_{n,m}^x+S_{i,j}^y{S}_{n,m}^y)]$$

$$⟶\sum _{i,j,n,m}{J}_{n,m}^{i,j}\left[⟨S_{i,j}^z⟩S_{n,m}^z+S_{i,j}^z⟨S_{n,m}^z⟩\right].$$

$$H_{\mathrm{MF}}$$

$$[H_{\mathrm{MF}},S^z]=0.$$

$$S_{2j+1}^a=-\frac{1}{2}{\tau }_{2j+1}^a,a=y,z,S_{2j+1}^x=\frac{1}{2}{\tau }_{2j+1}^x.$$

$$H_0$$

$$H^{\prime }$$

$$[H_{\mathrm{MF}},\sum _j(-1{)}^j{\tau }_j^z]=0.$$

$${\chi }^{\alpha \beta }(\omega ,Q)$$

$$S_{\mathrm{AFM}}^{xx}(\omega ,Q)$$

$$S_{\mathrm{AFM}}^{zz}(\omega ,Q)$$

$$|Q,n⟩=\frac{1}{\sqrt{L}}\sum _j{e}^{iQj}|j,n⟩.$$

$$⟨Q,m|H|Q,n⟩$$

$$G(j,k)=⟨Q,j|{\left(\omega -H_{\mathrm{MF}}+E_0+i\eta \right)}^{-1}|Q,k⟩.$$

$$|\mathrm{GS}⟩$$

$$⟨Q,n|{\tau }_Q^x|\mathrm{GS}⟩$$

$$⟨Q,n|{\tau }_Q^z|\mathrm{GS}⟩$$

$$⟨\mathrm{GS}|{\tau }_Q^z|\mathrm{GS}⟩$$

$${\chi }^{xx}(\omega ,Q)$$

$${\chi }^{yy}(\omega ,Q)$$

$${\chi }^{zz}(\omega ,Q)$$

Full text

I Introduction

Over the course of the last two decades, quasi one dimensional (Q1D) quantum magnets have been established as an ideal testing ground for key concepts of quantum many-particle physics such as quantum criticality QC ; QC2 ; QC3 , condensation of magnetic excitations affleck90 ; tsvelik90 ; affleck91 ; condensation1 ; condensation2 , quantum number fractionalization FT1 ; FT2 ; KCuF3 , dimensional crossover PhysRevLett.85.832 ; PhysRevB.71.134412 and confinement of elementary particles. Confinement originally arose in the context of high-energy physics as a pivotal property of quarks, but subsequently was realized to emerge quite naturally in one dimensional quantum many-particle systems and field theories featuring kink or soliton excitations Mccoy.PRD.18.1259 ; McCoyWu.PRB.18.4886 . The simplest such example involves domain wall (“kink”) excitations in Ising-like ferromagnets, and has been explored in exquisite detail in a series of experiments by Coldea and collaborators Coldea.Science.327.177 . Confinement in ladder materials was studied in Ref. Lake.Nat.phys.6.50, , while the confinement of spinon excitations has been recently investigated on the Q1D spin-1/2 Heisenberg-Ising antiferromagnetic compound BaCo2V2O8 Grenier.PRL.114.017201 . Here the spinon continuum, characteristic of 1D spin-chain, observed above the three dimensional ordering temperature $T_{N}$ that breaks up into a sequence of gapped, resolution limited modes in the 3D ordered phase ($T<T_{N}$). An interesting difference to the ferromagnetic case is that two sequences of bound states with longitudinal and transverse polarizations respectively have been observed.

In the present study we use inelastic neutron scattering to investigate magnetic excitations in the Q1D spin-1/2 XXZ system SrCo2V2O8 as a function of temperature covering both the 1D ($T>T_{N}$) and 3D ($T>T_{N}$) magnetic states. The experimental results are complemented by detailed theoretical considerations that provide a quantitative explanation of the experimental observations.

SrCo2V2O8 crystallizes in the centrosymmetric tetragonal space group $I4_{1}cd$ (No. 110) with lattice parameters $a=b=12.2710(1)$ Å and $c=8.4192(1)$ Å at room temperature Bera.PRB.89.094402 . The magnetic Co2+ ions are situated within CoO6 octahedra which form edge-sharing screw chains along the crystallographic $c$-axis Bera.PRB.89.094402 ; HePRB.73.212406 (Fig. 1(a)). There are four screw chains per unit cell which rotate in the $ab$-plane around (1/4, 1/4), (1/4, 3/4), (3/4, 1/4) and (3/4, 3/4) (Fig. 1(b)). Two diagonal chains rotate clockwise and the other two chains rotate anti-clockwise while propagating along the $c$-axis. This results in a complex interaction geometry with many possible superexchange interaction pathways. The strongest interaction is the antiferromagnetic intrachain coupling $J$ between nearest neighboring Co2+ ions along the chains. Weak interchain interactions are possible along both the sides ($J^{{}^{\prime}}$) and the diagonals ($J^{{}^{\prime\prime}}$) of the $ab$-plane (Fig. 1(b)). These interchain interactions are in fact probably comprised of several interactions due to the screw chain structure, some of which also have components along the $c$-axis as found for the isostructural compound SrNi2V2O8 Bera.PhysRevB.91.144414 .

The interchain interactions stabilize long-range collinear antiferromagnetic (AFM) order below $T_{N}=5.2$ K Bera.PRB.89.094402 with the spins pointing parallel to the $c$-axis (chain axis). Consecutive spins order antiferromagnetically along the chains while within the $ab$ plane, the spins order ferromagnetically/antiferromagnetically along the $a/b$ axis. The magnetic moment of the Co2+ ions in the distorted octahedral crystal field environment is described well by a highly anisotropic pseudospin, ${S}=1/2$ Abragram.book . The exchange interactions between the pseudospins in SrCo2V2O8 can be modeled by the Hamiltonian BonnerPR.135.A640

[TABLE]

where $S^{\alpha}_{i,j}$ is the alpha component of the $i^{th}$ spin of the $j^{th}$ chain. $J>0$ is AFM nearest neighbor intrachain exchange interaction and $J^{i,j}_{n,m}$ is the interchain interaction between the $i^{th}$ spin of the $j^{th}$ chain and the $n^{th}$ spin of the $m^{th}$ chain. The anisotropy parameter $0<\epsilon<1$, takes into account the XXZ-type anisotropy interpolating between the Heisenberg ($\epsilon=1$) and Ising ($\epsilon=0$) limits.

II Experimental Methods

Single crystals of SrCo2V2O8 were grown using the floating-zone method Bera.PRB.89.094402 . Inelastic neutron scattering (INS) experiments were performed using the cold neutron triple-axis-spectrometers FLEXX at Helmholtz-Zentrum Berlin, Germany, and PANDA at the Heinz Maier-Leibnitz Zentrum, Garching, Germany. Measurements were performed on a large cylindrical single crystal (weight $\sim 4.5$ g, diameter $\sim 4$ mm and length $\sim 40$ mm) in the ($h$,0,$l$) reciprocal space plane. The measurements were performed with fixed final wave vectors of $k_{f}$=1.3 Å*-1*, $k_{f}$=1.57 Å*-1* and $k_{f}$=1.8 Å*-1*. For these measurements, the sample was mounted on an aluminum sample holder and was cooled in cryostat. For the FLEXX spectrometer, a double focusing monochromator and a horizontally focusing analyzer were used. For the PANDA spectrometer, both monochromator and analyzer were double focusing. Higher order neutrons were filtered out by using a velocity selector on the FLEXX spectrometer and a cooled Beryllium filter on the PANDA spectrometer. Measurements took place at various temperatures between 0.8 K and 6.0 K.

III Experimental Results

III.1 High temperature phase $T>T_{N}$

Since for SrCo2V2O8 $T_{N}/J\ll 1$, we expect there to be a temperature regime $T_{N}\lesssim T\ll J$ in which the physics is essentially one-dimensional (1D) and approximately described by an anisotropic spin-1/2 Heisenberg XXZ chain. Single crystal inelastic neutron scattering measurements of SrCo2V2O8 at 6 K ($>T_{N}=5.2$ K) along the (0, 0, $l$)-direction (chain direction) reveal a gapped scattering continuum (Fig. 2a). For such wave vectors the polarization factors are such that only the components of the dynamical structure factor transverse to the direction of magnetic order contribute to the scattering cross section. The gap value of $\simeq 0.95$ meV at the (0,0,2) zone center, is quite small compared to the bandwidth of the dispersion ($\simeq 14.5$ meV) revealing that the compound lies intermediate between the Ising (gap $\sim J$, bandwidth $\simeq\epsilon J$) and Heisenberg (gapless, bandwidth $\simeq\pi J$) limits.

III.1.1 Spinon continuum at $T>T_{N}$

The observed spinon continuum at $T=6$ K ($\ll J$) is in good agreement with the predictions for the transverse dynamical structure factor of the integrable Heisenberg XXZ chain at zero temperature Caux2008 . The lower boundary of the two-spinon continuum as a function of reduced momentum transfer $0\leq Q<\pi$ along the chain is given by

[TABLE]

and the upper boundary is

[TABLE]

Here

[TABLE]

where $\kappa=\cos(Q_{\kappa})=\frac{1-k^{{}^{\prime}}}{1+k^{{}^{\prime}}}$, $I=JK(\sqrt{1-k^{\prime 2}})\sqrt{1-\epsilon^{2}}/\pi$, $K(k)$ is the complete elliptic integral of the first kind and the parameter $k^{\prime}$ is given by

[TABLE]

The specific value $Q_{\epsilon}$ is obtained from the solution of a quartic equation Caux2008 .

For SrCo2V2O8, $Q$ can be written in terms of the crystallographic wavevector transfer $Q^{\prime}$ as $Q=Q^{\prime}_{c}/4=2\pi l/4$, where $Q^{\prime}_{c}$ is the wave vector transfer in terms of the $c$-lattice parameter of SrCo2V2O8. The factor of four arises from the four equivalent Co2+ ions per unit cell along the chain direction ($c$-axis). Fitting the experimental continuum boundaries of SrCo2V2O8 to the above expressions yields the values of $J\approx 7.0\pm 0.2$ meV and $\epsilon\approx 0.56\pm 0.02$. The fitted continuum boundaries are represented by the solid black lines plotted over the data in Fig. 2(a).

Ref. Caux2008, also provides the theoretical expression for the transverse structure factor of the 2-spinon continuum. Using the fitted values of $J$ and $\epsilon$ for SrCo2V2O8, the calculated transverse structure factor is shown in Fig. 2(b) and can be directly compared to the experimental data in Fig. 2(a). Fig. 2(c) shows energy scans at several fixed wave vectors from (0,0,2) to (0,0,3) which pass through the lower edge of the continuum of SrCo2V2O8. The lines through the data are the theoretical intensities convolved with the instrumental resolution. Good agreement is achieved between experiment and theory except at (0,0,2) where the effects of interchain coupling and finite temperature which are not included in the calculation may alter the spectrum at lowest energies.

III.1.2 Villain mode

An interesting feature in the dynamical response of spin chains is the existence of a finite temperature resonance known as a Villain mode Villain . This “mode” was first observed by neutron scattering in Ref. Nagler1982, ; Nagler1983, and is a fairly general feature of spin chain models James.PhysRevB.78.094411 ; Goetze.Phys.Rev.B.82.104417 . The Villain resonance in the XXZ chain has been investigated theoretically by developing a perturbation theory around the Ising limit James.Phys.Rev.B.79.214408 . A prediction of this theory is that above a certain temperature, a narrow resonance develops at an energy

[TABLE]

The resonance corresponds to transitions between thermally occupied states and therefore disappears at zero temperature. In our case, we expect to see a resonance at low temperatures at

[TABLE]

which follows a similar dispersion to that of the lower boundary of the continuum but is shifted downward from it by an energy similar to the energy gap $\simeq$ 0.95 meV. The predicted Villain mode is indicated in Fig. 2(a) by the dashed black curve. As $T=6$K which is still quite low compared to the intrachain interaction ($T/J\approx 0.07$), we expect the temperature effects on the $T=0$ two-spinon continuum to be weak. Hence, the most noticeable effect of temperature is the emergence of additional peaks associated with the Villain mode in the $T=6$ K data just below the two-spinon continuum. A weak peak is indeed visible in the (0,0,3) data at $\approx 6.8$ meV and in the (0,0,2.75) and (0,0,2.5) scans at $\approx 6.3$ meV and $\approx 4.8$ meV, respectively (see Fig. 2(c)). These peak positions along with those obtained from other energy scans (not shown) are represented by the yellow circles in Fig. 2(a) and follow the predicted Villain mode dispersion given by the dashed black curve.

III.2 Low temperature phase $T<T_{N}$

Below its Néel temperature $T_{N}=5.2$ K, SrCo2V2O8 develops long-range magnetic order where the Co2+ spins order antiferromagnetically along the chains with their moments aligned parallel to the $c$-axis Bera.PRB.89.094402 . The dynamical structure factor well inside the ordered phase at $T=1.5$ K along the (1,0,$l$) direction is shown in Fig. 3. Its gross features including the total bandwidth and the energy gap are similar to those observed above $T_{N}$ (Fig. 2(a)). The weak scattering at (1,0,3) is associated with the fact that there are four equivalent screw chains per unit cell each with four Co2+ ions per $c$-lattice parameter. Neglecting interchain interactions this gives rise to a total of four “copies” of the cross section for a single chain, which are shifted with respect to one another by reciprocal lattice units along the chain direction (for details see Ref. Bera.PhysRevB.91.144414, on the isostructural compound SrNi2V2O8). For uncoupled chains we thus expect the intensity to be of the form

[TABLE]

As a result every reciprocal lattice point is an antiferromagnetic zone center for at least one of these copies, but their overall intensities $A_{l}({\bf Q}^{\prime})$ depend on the full momentum transfer ${\bf Q}^{\prime}$ and can be very different. For (1,0,$l$) all four independent contributions are present shifted consecutively by $\Delta Q_{c}^{\prime}=1$ r.l.u. along the chain. Their lower and upper boundaries are indicated by the different colored lines in Fig. 3(a). For (0,0,$l$) only a single contribution is visible, as observed in Fig. 2(a).

Careful inspection of the cross section at the antiferromagnetic zone center, reveals that the continuum observed at 6 K is transformed into a sequence of discrete, resolution-limited excitations at 1.5 K. As shown in Fig. 3(b) at wave vector transfer (1,0,2) nine peaks, labelled $E_{1}$-$E_{9}$, are observable in the energy range $\sim$ 1.5-5.5 meV. Since these discrete modes appear below the ordering temperature, they must arise from the interchain coupling. A detailed examination shows that each of the sharp peaks at (1,0,2) in fact consists of two closely spaced peaks with the higher energy peak being relatively weaker. For the wave vector (0,0,2) a single series of peaks is found (Fig. 3(c)), while at (3,0,1) both series of peaks are visible with similar intensities (Fig. 3(d)).

We have investigated the nature of the two series of peaks in more detail at several AFM zone centers with different wave vector components $Q_{a}^{\prime}$ and $Q_{c}^{\prime}$ (along the $a$ and $c$ axes) respectively. The measurements were performed over the lowest energy peaks around $\sim 1.5$ meV, see Fig. 4(a). The results indicate that when the wave vector transfer is parallel to the chain direction, e.g. (0,0,2), only one peak is present. If the $a$-component of the wave vector transfer is non-vanishing a second peak appears at higher energy. The relative intensity of the higher energy peak increases with increasing $Q_{a}^{\prime}$.

This intensity dependence provides important information about the nature of the two series of peaks. Neutron scattering is only sensitive to fluctuations perpendicular to the wave vector transfer. The higher energy series of modes that is absent for wave vector transfers parallel to the $c$-axis but becomes visible when $Q_{a}^{\prime}\neq 0$ must therefore be due to fluctuations along to the $c$-axis. We refer to this series of modes as longitudinal modes (L-modes) since they are due to fluctuations parallel to the ordered spin direction. In contrast, the available evidence suggests that the lower energy series of modes is associated with fluctuations in the $ab$-plane. We will therefore refer to these excitations as transverse modes (T-modes). The bound modes in SrCo2V2O8 were observed previously using terahertz spectroscopy as described in Ref. Zhe.PRB.91.140404, . This techniques allows the transverse modes to be measured to very high resolution, but the longitudinal modes are not accessible.

It is clear from Fig. 4 that the energies of the modes vary from one AFM zone center to another as a result of the interchain interactions. In order to investigate these interactions, the dispersions of the lowest energy pair of bound modes were measured along $Q^{\prime}_{a}$ by performing a series of energy scans at the constant-wave vectors ($h$, 0, 1) and ($h$, 0, 2) for various values of $h$ (Figs. 5(a) and 5(b)). Both the L-mode and T-mode disperse over a narrow bandwidth of $\leq 0.15$ meV, the modes are in-phase for $l$=2 but out-of-phase at $l$=1 (Fig. 5(c)). The dispersions are complex due to the many possible interchain interactions allowed by the screw-chain crystal structure (Fig. 1) that can reinforce or act against each other depending on the reciprocal lattice points as found for the isostructural compound SrNi2V2O8 Bera.PhysRevB.91.144414 . To fully quantify the strengths of the interchain interactions further measurements are required.

III.2.1 Modeling the energies of the observed bound modes

As we will detail in section IV, the bound states observed at low temperatures can be understood in terms of confinement of spinon pairs. The physical picture is that the interchain coupling induces a linearly confining potential between the elementary spinon excitations of the 1D chains. This was shown by Shiba in Ref. shiba, in the large anisotropy limit $\epsilon\ll 1$ of the model (1). In this limit spinons can be thought of as antiferromagnetic domain walls. As we will see in section IV the spinon confinement picture extends all the way up to the Heisenberg limit $\epsilon=1$. This suggests that the bound mode energies $E_{j}$ at the AFM zone center can be approximately extracted from the 1D Schrödinger equation describing the centre-of-mass motion of the spinon pairs

[TABLE]

Here $\mu$ is the reduced mass, $E_{0}$ is the spinon gap in absence of the confining potential, $\lambda$ is the molecular field at the Co2+ site produced by the interchain interactions, and the interaction potential between the two spinons is assumed to be a linear function of their separation $x$. The Schrödinger equation (8) has been previously applied successfully to describe aspects of confinement in the transverse field Ising chain Mccoy.PRD.18.1259 ; McCoyWu.PRB.18.4886 ; Isingconfinement1 ; BT ; Rutkevich2010 and in real materials Coldea.Science.327.177 ; Grenier.PRL.114.017201 . The solutions of Eq. (8) are given by Airy functions Landau and the corresponding bound state energies are

[TABLE]

where $\alpha=[\lambda^{2}(\hbar^{2}/\mu)]^{1/3}$ and the $\zeta_{j}$’s are the negative zeros of the Airy function. We use Eq. (9) as a phenomenological expression for the bound state energies, and fit the two parameters $E_{0}$ and $\alpha$ to our experimental data for the longitudinal and transverse modes separately. This gives excellent agreement with the observed spectra in all cases see Fig. 3(e) to (g). The fitted value of $\alpha$ is $\alpha=0.22\pm 0.1$ while the spinon gap $0.98<2E_{0}<1.21$ shows some variation between different AFM zone centers probably due to the interchain coupling.

It should be noted that the bound modes have their minimum at the reciprocal lattice points and disperse along the chain direction as can be observed in Fig. 3(a). The bound mode dispersion in the vicinity of the antiferromagnetic zone center is of the form

[TABLE]

where $Q$ is the reduced wave vector along the chain direction ($Q=Q^{\prime}_{c}/4$) and $m_{j}$ is the mass of the $j^{\rm th}$ bound state.

III.3 Temperature effects

In the simplest model the confining potential for spinons is proportional to the magnitude of the ordered moment and we therefore expect the bound modes to be sensitive to temperature at $T\approx T_{N}$. The temperature-dependence of the transverse and longitudinal bound spinon modes at the reciprocal lattice points (0,0,2) and (3,0,1) are shown in Fig. 6. As temperature approaches $T_{N}$ from below these modes broaden, become weaker and shift to lower energy. This shift is due to the weakening of the confining molecular field from the neighboring chains as the order moment value decreases with increasing temperature.

Another feature in the data is a strong broad peak centered around $E=0$ [Fig. 6(b)]. It is visible at $T\sim T_{N}$ but disappears for $T\ll T_{N}$ suggesting that it is due to short-range order between the chains that sharpens into magnetic Bragg peak position well below $T_{N}$. We observe that the peak is present at (3,0,1) but not at (0,0,2). The likely origin of this difference is that at (0,0,2) we only observe transverse correlations, while the peak is related to (emerging) 3D order along the longitudinal direction.

In addition to these changes, a sharp peak appears at (3,0,1) at the energy $E\sim 0.65$ meV for $T=4$ K and shifts towards lower energy with increasing temperature. No such temperature-induced peak is observed at (0,0,2), confirming that this feature is associated with the longitudinal structure factor. We attribute this peak to transitions between bound modes. At finite temperatures the lowest energy bound mode will become thermally populated and transitions between it and the higher energy bound modes are possible. Since the peak is observed in the longitudinal structure factor it arises from a transition between longitudinal bound modes or between transverse bound modes but not from a longitudinal to a transverse bound mode or vice versa. The distinction should be drawn between this feature which is a transition between thermally excited bound-spinon modes observed close to but below $T_{N}$ and the Villain mode which arises from transitions between thermally excited single spinon states which are observed above $T_{N}$.

If we denote the dispersion relation of the jth transverse/longitudinal bound state by $\Omega^{(a)}_{j}(Q)$ ($a=T,L$), transitions between them occur at energies

[TABLE]

In the low temperature ordered phase $\Omega^{(a)}_{j}(Q)$ will be $\pi$-periodic functions, so that the momentum transfer of the transition (11) will be $Q-P\ {\rm mod}\ \pi$. If we take $Q=P=\pi$ there will be transitions with momentum transfer zero and $\pi$. Setting aside the issue which transitions will give non-vanishing contributions to the dynamical structure factor (which is $2\pi$-periodic), we have checked whether the energy differences between bound modes at (3,0,1) are reflected in the energy of the temperature induced peaks. The inset of Fig. 6(b) shows the temperature-dependent energy of the thermally excited peak (green circles) which is in good agreement with the energy difference between the first and second longitudinal bound modes as a function of temperature (red squares). The above interpretation of temperature effects is supported by the theoretical analysis summarized in section IV.4.2.

IV Theory

As we have seen above, at $T>T_{N}$ the neutron scattering intensity is well described by a model of uncoupled spin-1/2 Heisenberg XXZ chains. At low temperatures interchain coupling effects are obviously important. In the following we constrain our analysis to a simple mean-field treatment of these interactionsshiba ; MFA1 ; MFA2 ; MFA3 ; review

[TABLE]

Using the fact that there is Néel order at low temperatures this leads to a description in terms of decoupled chains in a self-consistent staggered magnetic field

[TABLE]

The effective staggered field $h$ is a function of $J^{i,j}_{n,m}$ and temperature. We note that the Hamiltonian (13) has a U(1) symmetry of rotations around the z-axis

[TABLE]

In the following we analyze the dynamical structure factor in the model (13) by several different methods in various parameter regimes.

IV.1 Strong coupling expansion

A fairly comprehensive qualitative picture of the physical properties of the model (13) can be obtained by considering the strong anisotropy limit $\epsilon\ll 1$. This limit is amenable to an analysis by the method of Ishimura and Shiba IS and has been previously considered by Shiba shiba . As Ref. shiba, only considered the transverse component of the dynamical structure factor, we now give a self-contained discussion of this approach and then discuss the resulting picture for dynamical correlations. We find it convenient to map (13) to a ferromagnet by rotating the spin quantization axis on all odd sites around the x-axis by 180 degrees

[TABLE]

Here $\tau_{j}^{\alpha}$ are Pauli matrices. In terms of the new spins we have $H_{\rm MF}=H_{0}+H^{\prime}$ with

[TABLE]

where $\tau^{\pm}=\frac{\tau^{x}\pm i\tau^{y}}{2}$. The U(1) symmetry (14) gives rise to the commutation relations

[TABLE]

The zero temperature dynamical susceptibilities are given by

[TABLE]

where $\eta>0$ is infinitesimal, $E_{0}$ is the ground state energy and $\tau^{\alpha}_{Q}=\frac{1}{\sqrt{L}}\sum_{j}e^{iQj}\tau^{\alpha}_{j}.$ The dynamical structure factor of the antiferromagnetic spin chain (13) of interest is

[TABLE]

We will analyze (18) by carrying out a strong coupling expansion in the limit $\epsilon,h\ll J$ IS . Our starting point is the Ising part $H_{0}$ of the mean field Hamiltonian. The ground states of $H_{0}$ are simply the saturated ferromagnetic states $|\uparrow\rangle$ and $|\downarrow\rangle$ respectively. Their energies are $E^{(0)}_{0}=-\frac{JL}{4}.$ Spontaneous symmetry breaking selects e.g. $|\uparrow\rangle$. The low-lying excitations are then 2-domain wall states of the form $...\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\uparrow\uparrow\uparrow...$ We denote these by $|j,n\rangle$ where $j$ is the position of the first down spin and $j+n-1$ the position of the last down spin. The energies of these states are $E^{(0)}_{2}=-\frac{JL}{4}+J.$ A convenient orthonormal basis of 2-domain wall states with momentum $Q$ is obtained by taking appropriate linear combinations

[TABLE]

The matrix elements of the Hamiltonian in these states are

[TABLE]

Importantly, the only non-zero matrix elements in $H$ occur between domain-wall states with both even or both odd lengths. This is a consequence of the U(1) symmetry (17) and expresses the fact that acting with $H$ does not change the staggered magnetization (or equivalently the magnetization in the original $S^{\alpha}_{j}$ spin variables). Given (21) it is a straightforward matter to numerically compute the Green’s functions

[TABLE]

IV.1.1 Ground state in perturbation theory

The first order correction to the ground state is obtained by standard perturbation theory

[TABLE]

This gives the following matrix elements of spin operators

[TABLE]

IV.1.2 Dynamical Structure Factor

Substituting (24) into (18) then leads to the following approximate expression for the dynamical susceptibilities at $\omega>0$

[TABLE]

We note that these are consistent with the U(1) symmetry of rotations around the z-axis for the antiferromagnetic model (13). It is now straightforward to compute the dynamical structure factor (DSF) (19) numerically. Results for momentum transfer $\pi$ are shown in Figs. 7.

We see that the transverse DSF ($S^{xx}_{\rm AFM}(\omega,\pi)$) only ‘couples’ to half the bound states, while the longitudinal DSF ($S^{zz}_{\rm AFM}(\omega,\pi)$) is sensitive to the other half. This is in perfect correspondence with the experimental observations. The selection rule that gives rise to this behaviour is related to the conserved $S^{z}$ quantum number (14). It is clear from (21) that the Hamiltonian in the 2-domain wall sector is block diagonal in a basis of domain wall states of odd/even length. In terms of the original spins even/odd length domain walls correspond to even/odd values of the conserved $S^{z}$ quantum number (assuming the lattice length to be divisible by 4). This implies that there is one sequence of bound states with $S^{z}=0$, and a second with $S^{z}=\pm 1$. The first is visible in the longitudinal structure factor $S^{zz}_{\rm AFM}$, while the second contributes only to $S^{xx}_{\rm AFM}$. In the strong anisotropy limit we therefore have the simple cartoon picture for the physical nature of the bound states shown in Fig. 8.

IV.1.3 Gap as a function of field

The position of the first peak at $Q=\pi$, $E^{(T)}_{1}$, gives the excitation gap. Based on the relation of our problem to a Schrödinger equation with linear potential we expect

[TABLE]

In Fig. 9 we show the results obtained in our strong coupling expansion and a fit to (26), which is seen to give a very good account of the data [see Fig. 10].

IV.2 Field Theory in the vicinity of the isotropic point $\epsilon=1$

The physical picture obtained in the large anisotropy limit $\epsilon\ll 1$ remains valid in the entire regime $0<\epsilon<1$. To see this we consider the limit of weak anisotropy $\epsilon\approx 1$, where the mean-field Hamiltonian (13) can be written in the form

[TABLE]

where $\delta\approx J(1-\epsilon)$. In the parameter regime $h\ll\delta\ll J$ this model can be bosonized following e.g. boso ; boso2 , which leads to a two-frequency sine-Gordon model

[TABLE]

where $\lambda\propto\delta$ and $\mu\propto h$. The model (28) in the relevant parameter regime has been studied previously by a number of authors affleck ; DM ; BT . A fruitful line of attack is to start with ${\cal H}$ for $\mu=0$, and consider the $\mu$-term as a perturbation. The minima of the potential for $\mu=0$ and $\lambda>0$ occur at

[TABLE]

The solutions to the classical equations of motion are solitons and anti-solitons. Solitons interpolate between neighbouring vacua, e.g.

[TABLE]

while antisolitons have the opposite asymptotics

[TABLE]

At the quantum level solitons and antisolitons turn into elementary excitations of the sine-Gordon model.

IV.2.1 Soliton-antisoliton states

Following affleck ; DM we start with the soliton-antisoliton sector. We take the positions of the soliton and antisoliton to be $x_{1}$ and $x_{2}$ respectively and denote the classical energy for $\mu=0$ by $2\Delta_{s}$. When $\mu>0$, soliton-antisoliton states acquire an extra contribution to the energy

[TABLE]

In a non-relativistic approximation we then obtain a single-particle Schrödinger equation for the relative motion ($x=x_{2}-x_{1}$) with Hamiltonian

[TABLE]

Here $M_{s}=\Delta_{s}/v^{2}$ and the reduced mass is $\frac{M_{s}M_{\bar{s}}}{M_{s}+M_{\bar{s}}}=\frac{M_{s}}{2}$. This Schrödinger equation can be solved exactly in terms of Airy functions Landau , and the corresponding eigenstates describe the confinement on solitons and antisolitons. The bound state energies follow from the boundary conditions imposed on the wave function. If we require the wave function to be antisymmetric and therefore vanish at zero, we obtain

[TABLE]

Symmetric wave functions would instead lead to a spectrum of the form

[TABLE]

As soliton-antisoliton states have the same $S^{z}$ value as the ground state, the bound states (34) will be visible in the longitudinal structure factor.

IV.2.2 Soliton-soliton states

The considerations for two-soliton states are analogous. Classically the parameter $\mu$ characterizing the confining potential (32) is the same as in the soliton-antisoliton sector, but we don’t expect this to be true at the quantum level. We account for this by a different strength $\bar{\mu}$ of the potential, which then gives a sequence of energies

[TABLE]

Here we have taken the wave function to be antisymmetric because the zero-momentum limit of the soliton-soliton scattering matrix is $-1$. As soliton-soliton states have $S^{z}=1$, the bound states (36) will be visible in the transverse structure factor.

IV.2.3 Dynamical Structure Factor at $Q\approx\pi/a_{0}$

In the field theory limit the staggered magnetizations are given by

[TABLE]

Close to the antiferromagnetic wave number $\pi/a_{0}$ (where $a_{0}$ is the lattice spacing) the components of the DSF are thus given by

[TABLE]

In the longitudinal structure factor we therefore see confined soliton-antisoliton states, while the transverse components are sensitive to confined soliton-soliton and antisoliton-antisoliton states. Following Appendix B of DM we can derive expressions for the bound state contributions to the various correlators in leading order in perturbation theory in $\mu$ in the limit of very weak confinement

[TABLE]

where

[TABLE]

The function $f(\theta)$ is related to a particular two-particle form factor

[TABLE]

A similar analysis can in principle also be carried out on the level of the spin chain itself. The strengths of the confining potentials for spinon-antispinon and spinon-spinon two-particle states should be extracted from the known 4-particle form factor.

IV.3 DMRG Results

While the strong coupling expansion and field theory analysis provide a good qualitative picture of the dynamics, they do not apply quantitatively to the experimentally relevant regime $\epsilon\approx 0.56$. In order to overcome this shortcoming we have carried out numerical DMRG White.PRL.69.2863 ; Schollwock.AnnPhys.326.96 ; Hubig.PRB.91.155115 calculations with the SyTen tensor toolkit, based on the Hamiltonian (13). We first determine the gap of the lowest bound state in the $S^{z}=1$ sector as the difference between the lowest energies of the $S^{z}=0$ and $S^{z}=1$ sectors. This computation is extremely stable, even near criticality at small values of the field $h$ and with periodic boundary conditions. Finite-size effects can also be entirely removed by choosing sufficiently large system sizes. We are not able to determine the gaps of higher bound states in this way as this would require working at a fixed momentum. Fig. 10 gives the resulting values for the gap and a fit to the small-$h$ prediction (26). We see that already for system sizes $L=32$ the gap value is essentially converged and is in excellent agreement with the theoretical predictions for the scaling (26).

Second, we can calculate the ground state on a long chain, apply an excitation in the middle of the chain and then use Matrix Product State-based Krylov time evolutionmcculloch-krylov with matrix re-orthogonalization to evaluate the dynamical structure factors in the time-space domain. The Fourier transformation into momentum space is unproblematic. However, we are only able to evolve up to a time $t_{\mathrm{max}}J\approx 80$ for $\epsilon=0.2$. This limitation is a consequence of the entanglement growth during time evolution and the subsequent exponential increase in computational effort. This limit necessitates an articial damping factor $\mathrm{exp}\left(-\eta\frac{t}{t_{\mathrm{max}}}\right)$ to be introduced during the Fourier transform into frequency space. For $\epsilon=0.2$, $\eta\approx 1$ suffices and it is already possible to distinguish the physical peaks from the spectral leakage introduced by the transformation. For $\epsilon=0.5$, only slightly shorter time-scales $t_{\mathrm{max}}J\approx 50$ are achievable. Sufficient damping to remove spectral leakage then also removes the signal. To circumvent this problem, we use numerical extrapolation prior to the Fourier transformation to extend the data in time to very large $t^{\prime}_{\mathrm{max}}J\approx 1000$. We can then introduce a very small damping $\eta=1/700$ during the Fourier transformation and still remove all spectral leakage.

In Figs 11 we show results for the dynamical structure factor at momentum $Q=\pi$ on a system of $L=128$ sites with $h=0.05J$ and two values of the anisotropy $\epsilon$. In both cases we find two sequences of bound states associated with the transverse and longitudinal correlations respectively. The positions of the first peak in the transverse sector is in good agreement with Fig. 10.

IV.4 Tangent-space MPS methods

As we have explained in the previous section, targeting higher bound states variationally requires the ability to work within a fixed momentum sector. This is made possible by using tangent-space methods Haegeman2013b for matrix product states (MPS) that work directly in the thermodynamic limit. In particular, starting from a translationally-invariant MPS ground state on an infinite chain, we can apply the MPS quasiparticle ansatz Haegeman2012a to target the elementary excitations corresponding to isolated branches in the spectrum. This ansatz can be read as the MPS version of the Feynman-Bijl ansatz and single-mode approximation, but improves on these approaches in using the virtual degrees of freedom of the MPS to build an excitation on top of the ground state. As the ansatz explicitly contains a fixed momentum, it allows to systematically capture the wave functions of all quasiparticle excitations – the ones that contribute a $\delta$ peak in the DSF – within a certain momentum sector. In order to capture continuous bands in the spectrum, multi-particle excitations should be considered Vanderstraeten2015a ; Vanderstraeten2016 . Since the quasiparticle ansatz yields accurate variational expressions for the wave functions of the excited states, we can compute the energies and spectral weights for all states contributing to the DSF.

As it works directly in the momentum-energy plane, and does not suffer from finite-size effects, this method has access to the dynamical structure factor with perfect resolution. The only source of error is the variational nature of the approach, but the approximation can be systematically improved by growing the bond dimension of the MPS ground state. As the ansatz effectively exploits the correlations in the ground state to build an excitation, it can treat generic strongly-correlated spin chains with isolated branches in the spectrum to very high precision. An assessment of the accuracy of the variational wave function for a given excited state is provided by the variance $\langle H^{2}\rangle-\langle H\rangle^{2}$ of its energy, which can be evaluated exactly Vanderstraeten2015a .

At small values of the magnetic field $h$, we have a large number of stable bound states that live on a strongly-correlated background. Whereas time-domain approaches are necessarily limited in resolving the different modes, the quasiparticle ansatz is ideally suited for capturing all stable bound states with perfect resolution. Targeting the unstable bound states in the continuous bands would require a multi-particle ansatz Vanderstraeten2016 , but this has not proven to be necessary here. In order to compare with the experimental data, we have determined the dynamical structure factor by this method for several values of the anisotropy $\epsilon$ and the staggered magnetic field $h$. An additional shift of the energies has been introduced, corresponding to the three-dimensional dispersion of the modes. The best agreement with the experimental data is found for $\epsilon=0.56$ and $h=0.00643J$, and the corresponding structure factors are shown in Fig. 12.

IV.4.1 Comparison with strong-coupling expansion and DMRG results

It is useful to compare the results obtained by our different methods. We first consider a fairly strong anisotropy $\epsilon=0.1$ and weak field $h=0.025J$. Results for the quasiparticle ansatz (solid line) and the strong coupling approach (dashed line) are shown in Fig. 13. The agreement is seen to be good and any discrepancies can be attributed to the leading ${\cal O}(\epsilon^{2})$ corrections to the strong coupling result.

We have also compared the results of the quasi-particle ansatz to DMRG for $\epsilon=0.5$, cf. Fig. 14. The two methods are seen to agree very well. The remaining differences arise from the fact that the tangent-space MPS approach has been restricted to the calculation of the five lowest energy bound modes (in principle higher bound modes could be analyzed as well).

IV.4.2 Temperature effects

So far our theoretical analysis has been restricted to zero temperature. In order to access the $T>0$ regime we now combine the tangent-space MPS method with a low temperature linked cluster expansion of the dynamical susceptibility EK08 ; EK09 ; James.PhysRevB.78.094411 ; Goetze.Phys.Rev.B.82.104417 . The basic idea is to treat the low-temperature regime as a gas of bound states that scatter purely elastically. This is expected to be a good approximation as long as the temperature is small compared to the minimal gap $\Delta_{\rm min}$ of the lowest energy bound state, i.e.

[TABLE]

The main temperature effects are a broadening of the $T=0$ coherent single-particle peaks and the emergence of additional peaks in the dynamical structure factor, which correspond to transitions between thermally populated single-particle excitations. The first effect requires an analysis of matrix elements between single particle and two particle excitations. This is a non-trivial task beyond the scope of the present work. At sufficiently low frequencies $\omega<\Delta_{\rm min}$ the second effect is easier to capture. Let us denote the single-particle excitations of the $a^{\rm th}$ bound state with momentum $p$ by

[TABLE]

and the corresponding dispersion relations by $\epsilon_{a}(p)$. Then the leading contributions to the dynamical structure factor at low temperatures and frequencies are

[TABLE]

In Fig. 15 we show the contributions (44) due to transitions between thermally excited bound modes for the experimentally relevant parameter set $\epsilon=0.56$, $J=7{\rm meV}$, $h=0.00614J$ and $T=4K$. Transitions occurring at very low frequencies have not been taken into account, because the low energy regime is dominated by the broadened Bragg peak. For comparison the longitudinal (red dashed line) and transverse (solid blue line) components of the dynamical structure factor at $T=0$ are shown as well. A finite temperature resonance in the longitudinal structure factor at a frequency $\omega\approx 0.1J$ is clearly visible, while contributions to the transverse structure factor are very small.

V Theory vs experiment

We are now in a position to compare theoretical and experimental results. The first task is to determine appropriate parameters for applying the effective 1D model (13) to SrCo2V2O8. Estimates for the exchange $J$ and anisotropy $\epsilon$ were obtained in section III.1 by comparing the data collected for $T>T_{N}$ to the zero temperature transverse dynamical structure factor for (13) with $h=0$. Such a comparison is appropriate because $T\ll J$ and gives values of $J\approx 7.0$ meV and $\epsilon=0.56$. The remaining parameter is the strength $h$ of the effective staggered field. As this arises from a mean-field decoupling of the interchain interactions, it is temperature dependent. As shown in section IV, $h$ can be fixed by computing the energies of the first few bound states and comparing them to the measured peak positions. One caveat is that the gap of the lowest bound state is not necessarily well described by the purely 1D model (13). Indeed, in simple quasi-1D systems of weakly coupled chains corrections to the simple mean-field approximation due to the interchain couplings can be taken into account by a random-phase approximation, which gives the following expression for the dynamical susceptibility

[TABLE]

Here $J_{\rm int}({\bf q})$ is the Fourier transform of the interchain coupling, and we have assumed that we are dealing with a system of equivalent chains. It is clear from (45) that at a given wave vector the singularities of $\chi_{\rm 3D}(\omega,{\bf q})$ are shifted in energy by a constant compared to those of $\chi_{\rm 1D}(\omega,q)$. The situation in SrCo2V2O8 is much more complicated, because there are several counter rotating screw chains per unit cell. A refined analysis of the interchain coupling by a generalization of (45) is possibleZheludev but beyond the scope of this work.

Keeping this discussion in mind, we first try to obtain an optimal description of the energy splittings between the observed coherent modes by a pure one dimensional model. At temperature $T=1.5$ K we can reproduce the energy differences of the first few peaks with a value of $h\approx 0.00643J$. The resulting comparison between the transverse modes calculated by this mean-field model using the tangent-space MPS method (13) and the experimental data at $Q=(0,0,2)$ is shown in Fig. 16. We see that the mean-field model reproduces the experimental results very well up to an overall shift of about $0.122{\rm meV}$ in energy. Since the gap is very sensitive to corrections to the mean-field model, as can be seen from the RPA expression (45), such a shift is not surprising. Furthermore in the experimental data the interchain couplings give rise to a dispersion of the gap with a comparable bandwidth of $\sim 0.15$ meV (see Fig. 5).

The dynamical structure factor calculated in the mean-field model by the tangent-space MPS method was compared to the data at several reciprocal lattice points as shown in Fig. 17. Both the transverse and longitudinal structure factors are plotted and the effect of interchain coupling is taken into account by introducing a wavevector-dependent energy shift. At each wavevector the intensity of the two structure factors are weighted by their respective polarization factors due to the component of their magnetization perpendicular to the wavevector transfer (see section III.2) as well as by the square of their $g$-factors Zhe.PRB.94.125130 . An overall scaling factor is also introduced to match the theoretical intensity to the data and the theoretical peaks are convolved by a Gaussian to model the experimental resolution. The solid black line gives the sum of the two structure factors as well as a linear background and represents the expected neutron scattering intensities. Considering the highly complex counter-rotating screw chain structure and the many possible interchain interactions that are neglected in this calculation, the agreement between experiment and theory is remarkably good.

The finite temperature results of section IV.4.2 are also in good agreement with the experimental observations [Fig. 18]. We saw that at a temperature of $4$ K our one dimensional model displays a finite temperature resonance in the longitudinal structure factor at an energy of about $0.7$ meV. This is in good agreement with the experimental observation of a resonance in the longitudinal structure factor at $\omega\approx 0.65$ meV [Fig. 6(b)].

VI Discussion

We have presented results of inelastic neutron scattering experiments on the quasi-one dimensional spin-1/2 Heisenberg magnet SrCo2V2O8. Above the Néel temperature $T_{N}\approx 5.2K$, the neutron scattering cross section is dominated by a scattering continuum that is well described by a spin-1/2 Heisenberg XXZ chain with antiferromagnetic exchange $J\approx 7.0\pm 0.2\ {\rm meV}$ and anisotropy parameter $\epsilon\approx 0.56\pm 0.02$. The scattering continuum is formed by fractionalized $S^{z}=\frac{1}{2}$ spinon excitations. At temperatures below $T_{N}$ the structure factor exhibits two sequences of resolution-limited dispersing peaks that are associated with fluctuations along (L-modes) and perpendicular (T-modes) to the ordered magnetic moment respectively.

The origin of these coherent modes can be understood by a one dimensional model [Eq. (13)], in which a (temperature dependent) staggered magnetic field is generated in the ordered phase through a mean-field decoupling of the interchain interactions. The model can be studied analytically for both strong and weak exchange anisotropies and in both limits the effect of the staggered field is to confine the spinon excitations into two sequences of bound states. At intermediate values of the exchange anisotropy we have used DMRG and MPS methods to obtain quantitative results for the dynamical structure factor. It turns out that the experimentally relevant parameter regime cannot be reached even by state-of-the-art DMRG methods. Due to entanglement growth, the time scale by which dynamical correlation functions that can be computed by DMRG is restricted, which in turn imposes limitations on the achievable energy resolution. We therefore have employed a recently developed tangent-space MPS method, which is based on constructing MPS representations for excited states. Application of this method allows the computation of the dynamical structure factor, which is found to be in good agreement with experiment.

Our work establishes SrCo2V2O8 as a beautiful paradigm for spinon confinement in a quasi-one dimensional quantum magnet. There are a number of interesting questions that deserve further investigation. On the theoretical side a more involved investigation of the dynamical structure factor at finite temperatures would improve our understanding of the thermally induced peaks observed in the data both below $T_{N}$ (transitions between bound modes) and above $T_{N}$ (the villain mode). On the experimental side, the precise form of the interchain interactions needs to be clarified by extensive measurements of the bound mode dispersion relations perpendicular to the chain direction at lowest temperatures. We have seen that it is necessary to account for these interations beyond a simple mean-field decoupling in order to describe the data. As the crystal structure is rather complex this goes beyond the scope of the present work. Finally, it would be interesting to analyze the effects of an applied uniform magnetic field. Terahertz spectroscopy measurements reveal the emergence of novel excitations as a function of both transverse and longitudinal magnetic field Zhe.PRB.94.125130 ; Zhe_Long_Field which could be investigated using a combination of neutron scattering and the theoretical methods described here. We hope to return to these questions in future work.

Acknowledgements.

We acknowledge the Helmholtz Gemeinschaft for funding via the Helmholtz Virtual Institute (Project No. HVI-521). This work was supported by the EPSRC under grant EP/N01930X/1 (FHLE), the ExQM graduate school and the Nanosystems Initiative Munich (CH). L. Vanderstaeten acknowledge the financial support from FWO Flanders.