Tuning the confinement potential between spinons in the Ising chain CoNb2O6 using longitudinal fields and quantitative determination of the microscopic Hamiltonian

The Ising chain realizes the fundamental paradigm of spin fractionalization, where locally flipping a spin creates two domain walls (spinons) that can separate apart at no energy cost. In a quasi-one-dimensional system, the mean-field effects of the weak three-dimensional couplings confine the spinons into a Zeeman ladder of two-spinon bound states. Here, we experimentally tune the confinement potential between spinons in the quasi-one-dimensional Ising ferromagnet CoNb2O6 by means of an applied magnetic field with a large component along the Ising direction. Using high-resolution single crystal inelastic neutron scattering, we directly observe how the spectrum evolves from the limit of very weak confinement at low field (with many closely-spaced bound states with energies scaling as the field strength to the power 2/3) to very strong confinement at high field (where it consists of a magnon and a dispersive two-magnon bound state, with a linear field dependence). At intermediate fields, we explore how the higher-order bound states disappear from the spectrum as they move to higher energies and overlap with the two-particle continuum. By performing a global fit to the observed spectrum in zero field and high field applied along two orthogonal directions, combined with a quantitative parameterization of the interchain couplings, we propose a refined single chain and interchain Hamiltonian that quantitatively reproduces all observed dispersions and their field dependence.


I. INTRODUCTION
Fractionalization of coherently propagating spin-flips into two or more quasiparticles is a phenomenon of much interest in condensed matter physics.While in two and higher dimensions such phenomena require highly frustrated interactions, in one dimension, reduced mean field effects lead to fractionalization even in unfrustrated systems.A canonical example is the Ising chain in tilted field, with a deceptively simple Hamiltonian with J > 0 the Ising exchange, and h x and h z applied transverse and longitudinal fields, respectively.Consider first the case h z = 0: starting from a ferromagnetic alignment of all spins along the Ising z axis and flipping a single spin creates two domain walls.These can separate at no energy cost and move apart independently of each other in the presence of the transverse field [1,2], as illustrated in Fig. 1A.Therefore, a local spin flip, created for example in a neutron scattering process, fractionalizes into a pair of domain wall quasiparticles (spinons, also referred to elsewhere as solitons or kinks).However, in the presence of a finite longitudinal field h z > 0, domain walls are no longer free but have an attractive interaction, because there is an energy cost proportional to h z that increases linearly with their separation.The longitudinal field acts as an effective string tension between the domain walls, not allowing them to separate but confining them into bound states, realizing a one-dimensional analogue of the confinement of quarks into mesons [3].
By tuning the strength of the longitudinal field h z , one could then explore the crossover from the regime of weak confinement (with h z a small perturbation on the scale of J) where many closely-spaced bound states are expected, with energy separation predicted to scale as a power-law h 2/3 z , to the regime of strong confinement (h z comparable to J), where, depending on the relative sizes of J, h z and h x and on what other subleading exchange terms may be present in the spin Hamiltonian beyond the minimal model in (1), only one or at most two bound states are expected, with their energies scaling linearly with h z .The motivation behind the present studies was to explore the manifestation of this physics experimentally in a material where an external magnetic field can be used to tune the longitudinal field h z to cover the full range from weak to strong confinement.
The material CoNb 2 O 6 is considered to be an excellent experimental realization of a ferromagnetic Ising chain with a low enough exchange energy scale that the full phase diagram in magnetic field is experimentally accessible [4][5][6][7][8][9][10].It displays a quantum phase transition in transverse field, from an ordered phase to a quantum paramagnet, and around the critical point it displays the expected universal properties, such as evidence for an emergent E8 spectrum [4,7].The magnetic ions are Co 2+ , arranged into zigzag chains running along the crystallographic c-direction, with the zigzag in the bdirection, as shown in Fig. 2. The chains form a distorted triangular lattice in the ab-plane and the crystal structure is orthorhombic (space group P bcn) with lattice parameters a = 14.1337Å, b = 5.7019 Å and c = 5.0382 Å at 2.5 K [11].The combination of crystal field effects and spin-orbit coupling leads to an effective S = 1/2 Kramers In a finite longitudinal field hz, there is an energy cost linear in the separation, as if the two domain walls interacted via a string tension λ (curly line between solid dots), which stabilizes confinement bound states.C) In the weak confinement regime (hz ≪ J), many bound states exist, which are coherent superpositions of states with the two domain walls separated by many sites.Solid lines are the Airy wavefunctions for the first three bound states as per (4) for parameters λµ/ℏ 2 = 0.072, relevant for the longitudinal mean field due to interchain interactions in CoNb 2 O 6 [4].The domain wall separation, x, is in units of the spacing along the chain c/2.doublet ground state with strong Ising-like character, which is separated from the next lowest Kramers doublet by 30 meV [12].Weak interactions between chains stabilize magnetic order at low temperatures, and below 1.97 K, the spins are ferromagnetically aligned along each zigzag chain due to a dominant nearest-neighbour Ising exchange, with an antiferromagnetic pattern between chains [11,13].The magnetic moments lie in the ac-plane [14,15], at an angle γ = 30 • to the c-axis [11], which we take to be the local Ising z-axis; a field applied along the b-axis is then transverse to the Ising axes of all spins.The ferromagnetic nature of the dominant interactions makes CoNb 2 O 6 an ideal candidate for experimental tuning of the confinement potential between spinons.
Here, we report high resolution single crystal inelastic neutron scattering (INS) measurements that probe the full wavevector-dependence of the spectrum as a func-tion of magnetic field applied along the crystallographic a-axis, with a large longitudinal (along Ising axis) component.Our results complement earlier terahertz (THz) spectroscopy measurements for the same applied field direction, which probed the spectrum at the zone centre (wavevector transfer Q = 0) [9,16], as well as THz spectroscopy and INS measurements previously taken on the Ising material CoCl 2 • 2 H 2 O where the interchain couplings are more significant [17,18].In particular, we parameterize the evolution of the bound state spectrum with field and find good agreement with scaling laws expected in the regimes of weak and strong confinement, corresponding to low and high applied fields, respectively.We also parameterize the effects of the interchain couplings on the high-field dispersive modes for fields along both a and b to extract pure one-dimensional (1D) dispersion relations.We then compare those with exact diagonalization (ED) calculations to arrive at a refined spin Hamiltonian that quantitatively reproduces all the features seen in the full wavevector dependence of the inelastic neutron scattering data, across a wide range of applied fields.
The rest of this paper is organized as follows.Section II outlines the inelastic neutron scattering experiments and Sec.III presents the experimental results for the spectrum in field applied along a across a wide range of field strengths from weak confinement (III A) to strong confinement (III B), with the evolution of the high-order bound states presented in Sec.III C. Next, in Sec.IV, we propose a Hamiltonian for the in-chain interactions and refine the parameters by comparison to the observed dispersions at zero field and high field along two orthogonal directions (a and b).Sec.V contains our conclusions.The Appendices contain the characterization and parameterization of the effects of the interchain interactions, as well as comparisons between the data and other parameter sets.

II. EXPERIMENTAL DETAILS
Inelastic neutron scattering measurements of the magnetic excitation spectrum were performed using the direct geometry time-of-flight spectrometer LET at the ISIS facility [19].The sample was a large single crystal (4.59 g) of CoNb 2 O 6 grown by the floating-zone technique, as described in [20], and mounted in the (0kl) horizontal scattering plane, where the wavevector transfer Q is labelled as (hkl) in reciprocal lattice units of the structural orthorhombic unit cell, so Q = 2π h a , k b , l c .The sample was mechanically fixed in place using copper brackets so that it would not move or rotate due to the torques from the applied field.The sample was cooled using a dilution refrigerator insert and all data were collected below 0.14 K.A magnetic field up to 9 T was applied vertically, along the crystallographic a-axis, which has a longitudinal field component h z = g z µ B B sin γ, where B is the externally applied field magnitude.
LET was operated to measure simultaneously the inelastic scattering of incident neutrons with energies of E i = 2.46, 4.30 and 9.33 meV; the measured energy resolutions (full width half maximum, FWHM) on the elastic line were 0.038(1), 0.085(1), and 0.274(2) meV, respectively.To obtain the overview of the field dependence of the spectrum, the scattering was measured at each field for one fixed sample orientation chosen such that data projected along the chain direction covered a large part of the along-chain Brillouin zone.For these measurements, the sample was aligned with the (010) axis at an angle of 8 • to the incident beam direction, and scattering was measured for a typical counting time of 3 hours per field at an average beam current of 40 µA of protons on target.In addition to the single orientation measurements, multi-angle (Horace) scans were performed to obtain a full four-dimensional data set of scattering intensity as a function of three momentum directions and energy at two representative fields, 1.5 and 8 T. For these, the sample was rotated about the vertical a axis in steps of 2 • over an angular range of 88 and 108 • , with 11 and 17 minutes counting per step, respectively.Additional Horace scans were collected in a separate experiment on LET at 0.14 K in a field of 9 T ∥ b on the same crystal as in [21] mounted in the (h0l) horizontal scattering plane.For these measurements, the incident energies used were E i = 2.14, 4.02 and 10.17 meV, and the angular range was 145 • covered in steps of 1 • with 6 minutes counting per step.The Horace scans were used to quantitatively parameterize the interchain dispersions, as detailed in Appendix A. The raw time-of-flight neutron data were converted to scattering intensities S(Q, ω) using mantid [22], and the resulting data were then analysed using the horace [23] and mslice [24] packages.

A. Weak confinement regime
In the limit of small longitudinal field h z ≪ J, the confinement of domain walls into bound states can be captured via a Schrödinger equation for the domain wall separation distance x in the continuum limit [3,25,26] where µ is the reduced mass for the domain wall pair, λ = 2h z /c is an effective string tension between the domain walls, c/2 is the spacing of sites along the chain, and m 0 is the energy cost to create a single domain wall.For the Hamiltonian (1), m 0 = J/2 − h x [27].Eq. ( 2) can be recast as the Airy equation, giving energy eigenvalues where −z n are the zeros of the Airy Ai function; the corresponding eigenstates are The wavefunctions for the lowest three bound states are sketched in Fig. 1C.This shows that the bound states in the weak confinement regime are expected to be superpositions of many states, with finite amplitude even for states with large separations between domain walls.Eq. ( 2) is applicable for CoNb 2 O 6 even in the absence of an externally applied field.This is because i) there is a finite longitudinal (internal) mean field due to the weak interchain interactions h z = 2Jλ MF ⟨S z ⟩, with J and λ MF defined in Sec.IV A and ⟨S z ⟩ the expectation value of the spin component along the Ising axis; and ii) there are terms in the spin exchange Hamiltonian beyond the dominant Ising exchange (as will be shown in Sec.IV) that lead to domain wall propagation, so a finite µ.The bound state spectra observed in previous INS [4] and THz spectroscopy experiments [7,9] have shown good agreement with the predictions of (2) with h z interpreted as the interchain longitudinal mean field, which has the same magnitude for all sites in the zerofield antiferromagnetic phase.Zero-field data collected in the current experimental setting are shown in Fig. 3A, where the three lowest bound states are clearly visible at the lowest energies.Also notable is the sharp mode at the top of the spectrum with a dispersion curving the opposite way to the lower modes.This mode, which is stabilized by a different mechanism from the low-energy confinement bound states, is a single-spin-flip bound state stabilized by a subleading term in the spin Hamiltonian of the form −Jλ S j S + j S − j+1 + S − j S + j+1 /2.This XY exchange term allows a single-spin-flip state (i.e., two domain walls on adjacent sites) to hop between nearestneighbour sites along the chain as a single entity and this state was therefore dubbed a kinetic bound state [4].
The zigzag nature of the magnetic chains illustrated in Fig. 2 leads to a doubling of the unit cell along the chain direction compared to straight chains, therefore a doubling of the number of magnetic modes, so each "primary" mode with dispersion ω(h, k, l) has a "shadow" version with dispersion ω(h, k, l + 1) obtained from the "primary" version by Brillouin zone folding.In zero field, in the approximation of decoupled chains, the intensities of these modes are proportional to cos 2 (2πkζ) and sin 2 (2πkζ) respectively, where the magnetic ions alternate in position along the b direction as ±ζb between consecutive sites along the zigzag chain (ζ = 0.165) [21].In Fig. 3A, the top mode (the kinetic bound state) is a shadow mode with finite intensity only because of finite k, whereas all low-energy confinement modes are primary.Throughout this paper, we refer to the nth lowest energy (primary) mode as m n , in both the weak confinement and strong confinement regimes.Modes with the same label in these two regimes have very different character, but there is a smooth crossover between the two regimes.The remaining panels in Fig. 3 show the evolution of the spectrum upon increasing field.It is known that low fields applied along a induce a series of spin-flip transitions such that above a threshold field of 0.14 T, all spin components along a are parallel to the applied field [28][29][30].All finite field INS data were therefore collected at fields at and above this threshold field to ensure all spin sites experience the same magnitude longitudinal and transverse fields.Upon increasing applied field, all bound states move up in energy, with the higher order ones moving at a faster rate.This is because the higher order bound states contain in their superposition more weight for states with domain walls further apart [see Fig. 1C], so with more spins flipped opposite to the applied field direction, and therefore higher Zeeman energy or higher effective g-factor.Since the higher modes increase in energy more quickly, the relative energy separation between adjacent modes increases and modes become better resolved, such that whilst at 0.14 T (panel B) only m 1 to m 3 are resolved, at 0.5 T (panel D), m 1 to m 7 can be resolved.The kinetic bound state, however, increases in energy relatively slowly since it is a single-spin-flip state, such that the confinement bound states move past it as field is increased.Upon increasing field, the extent of the spectrum covers a progressively wider energy range and so, to probe the full spectrum above 2 T, we use a higher incident energy, coarser-resolution configuration with data plotted in Fig.

B. Strong confinement regime
The INS spectrum at 8 T, representative of the strong confinement regime, is shown in Fig. 5A.Two dispersive modes, lower energy m 1 and higher energy m 2 , are seen, each with both a (stronger intensity) primary and a (weaker intensity) shadow version (obtained from the primary via an l → l + 1 translation).In this regime, m 1 is a magnon mode, a coherently propagating single spin flip, as expected for the lowest energy excitation in a field-polarized paramagnet.Meanwhile, the m 2 mode can be understood as a dispersive two-spin-flip state, i.e. a bound state of two m 1 magnons.To understand why a two-magnon state can still exist in this limit of high fields, and why it is dispersive, albeit with a much suppressed bandwidth compared to m 1 , it is insightful to consider a minimal Ising-like XXZ model in longitudinal field The spectrum can be solved exactly [31] and contains, in addition to conventional one-magnon (m 1 ) excitations,  3), collected using a high-coverage coarser-resolution configuration (Ei = 9.33 meV) for the same fixed sample orientation as in Fig. 3, with wavevectors projected along the chain direction.Panel A is at the same field as the higher-resolution data in Fig. 3H.also a two-magnon (m 2 ) bound state, with dispersion relations and (un-normalized) wavefunctions, expressed in wavevector units relevant here, as and For m 1 , the dispersion comes at first order in the XY exchange Jλ S .For m 2 , the hopping process occurs in two stages via an intermediate state in which the spin flips are separated by one site [second term in (8)], therefore the hopping process is second order in the XY exchange, which leads to the dispersive term proportional to λ 2 S in (6).We propose that the m 2 state seen experimentally has the same qualitative content as the two-spinflip bound state for the minimal model in ( 5), but with modifications appropriate for the relevant Hamiltonian, such that the m 2 dispersion relation depends on all λ's in (9), to be discussed later.The dynamical correlations calculated for this model are shown in Fig. 5B and capture all key features of the dispersions and intensities of both m 1 and m 2 modes.
Note that a two-spin-flip state cannot usually be seen in inelastic neutron scattering [it would not be observable for the Hamiltonian in (5)], since the scattering intensity of a given state is determined by its overlap with a single spin flip created in a neutron scattering process.However, in the present case, the external field is applied along the crystallographic a-axis, which makes a finite angle 90 • − γ with the Ising direction, and so the field is not perfectly longitudinal but also has a finite transverse component.In addition, as we discuss in Sec.IV, there are also sub-leading exchange terms in the spin Hamiltonian that break spin conservation.Both the finite (large) transverse field and the (relatively much smaller) exchange terms in the Hamiltonian that break spin conservation lead to some mixing between states with different numbers of spin flips, allowing the predominantly two-spin-flip m 2 state to have a finite admixture of a single-spin-flip state and therefore to be visible in the present INS experiments.
It is notable that the m 2 mode is expected to survive up to indefinitely large field (although the intensity would become progressively weaker as the mixing with the single-spin-flip state would progressively decrease upon increasing Zeeman energy).This is because the binding energy of m 2 relative to the 2m 1 continuum for the minimal model in (5,6) is m 2 − 2m 1 = −J in the limit λ S → 0, i.e. the two magnons gain Ising energy if they are next to one another.The fact that the m 2 mode is lower in energy than the 2m 1 continuum means that it is not possible for it to decay, hence it survives as a sharp mode.
It is also important to note that at these high fields the m 1 state has non-negligible dispersion bandwidth perpendicular to the chain direction because of the finite interchain couplings; this is discussed further in Appendix A. This has the consequence that the lineshape of the m 1 mode in Fig. 5A appears artificially broadened because the data are integrated over a large range in the wavevector direction transverse to the chains.The interchain dispersion at high field is non-negligible because the hopping strength between chains for a single spin flip is first order in the interchain exchange.In contrast, at lower fields, even the lowest energy mode (m 1 ) has multiple spin flips, as illustrated in Fig. 1C for ψ 1 , so the interchain dispersion is higher order and is thus much suppressed, such that it is essentially undetectable at low field.
We also note that the bottom of the m 1 dispersion is visibly flattened compared to a perfect sinusoid [see Fig. 5A], indicating a double Fourier component along l, suggesting hopping to next-nearest-neighbour sites along the chain, which will be further discussed in Sec.IV A.

C. Fate of higher-order bound states
The evolution of the INS spectrum upon increasing field from 2 to 9 T is shown in Fig. 4. A notable feature is that, upon increasing field, the higher order bound states progressively become less dispersive and rapidly lose spectral weight, with the highest ones doing so at the fastest rate.In detail, at 2 T (panel A), m 1 and m 2 have similar intensities and bandwidths, m 3 is fairly flat and m 4 is very flat and fairly weak.Upon increasing field, m 4 moves up in energy and becomes weaker until it disappears between 4 and 4.5 T. At those fields, m 3 is completely flat and weaker in intensity, then becomes undetectable above 6.5 T.
The field dependence of the m 3 mode is summarized in Fig. 6A.The intensity (red triangles and right vertical axis) drops off very rapidly upon increasing field, decreasing by more than an order of magnitude between 1 and 4 T.In this field range, the m 3 energy (blue round symbols, left axis) is below the estimated 2m 1 continuum lower boundary (dashed line), so the m 3 mode should be stable, with no mechanism for decay; nevertheless, its intensity drops off very fast upon increasing field.A qualitatively similar behaviour is observed for the m 4 mode (panel B): the intensity drop-off is an even steeper function of increasing field and again the intensity decreases very quickly even before 2 T, above which m 4 enters the 2m 1 continuum.Within the sensitivity of the experiments, no systematic intensity change occurs for either of the m 3,4 modes as they enter the continuum and their signal could be followed even up to significantly higher fields, when the modes should be deep into the continuum Intensities are raw neutron scattering intensities on an arbitrary scale (the same for all panels).Note that a clear signal at the expected m3,4,5 energy is present even after crossing the 2m1 continuum lower boundary (dashed vertical line in panels C-F), which has been calculated using fitted dispersion relations for the m1 mode and assuming that m1 particles do not interact with each other.
energy range [see Fig. 6C-F].This might suggest that the matrix elements for the decay processes m 3,4 → m 1 + m 1 are relatively weak, and that the dominant effect leading to the fast intensity drop-off with increasing field is the progressively reduced mixing of the single-spin-flip state into the wavefunction of the m 3,4 modes, as upon increasing field the number of spin flips becomes a progressively better approximation for a good quantum number.This picture is consistent with the results of ED calculations, which indicate that while there is some decay of the bound states when they overlap with the continuum, the predicted broadening is of order 0.01 meV, which is not resolvable within the experimental results presented here.
As well as studying the disappearance of the higher bound states, the crossover between the weak and the strong confinement regimes was quantitatively tested by doing fits to the expected field dependence of the mode energies.In the high field regime, a linear dependence of mode energy on field is expected.This is because, in this regime, the number of spin flips in the mode is approximately a good quantum number, with the m n mode containing states with n spins flipped compared to the field-polarized state.Thus, the Zeeman energy is expected to scale linearly with both mode number and field.Indeed, the band averages are well fit by a linear dependence of energy on field and mode number, i.e., a fit where the gradients of the m 2 , m 3 and m 4 lines are constrained to be 2, 3 and 4 times respectively the gradient of the m 1 line, as illustrated in Fig. 7A.A cutoff field of 3.5 T, assumed to be close enough to the high field limit, was used, and only data at fields above this were included in the fit, which was performed simultaneously to the four lowest energy modes.The quality of this fit, especially for the m 1 and m 2 modes, also adds evidence that these are single-spin-flip and two-spin-flip modes respectively, i.e. derived from modes with wavefunctions given in (7) and (8).The fit in this regime used the band averages, as we have found that the band minimum in this regime is affected by the tilting of the local magnetization towards the field, which leads to a field-dependent bandwidth.This effect is quantitatively understood and A B FIG. 7. A) Band averages fit (solid lines) to a form En = ngµBB + E0n where n labels modes, g = 4.20( 2) is an effective g-factor (for fields along a) and E0n is an effective zerofield energy.The fit only includes data for B ≥ 3.5 T. B) One-dimensional band minima, corrected for the effects of interchain hopping, as a function of field in the low field regime.The band minima have been simultaneously fit (solid lines) to a form En = αzn|B − B0| 2/3 + E0 as per (2), where −zn are the zeroes of the Airy Ai function, B0 = −0.24(2)T is a field offset due to the mean field effects of the antiferromagnetic interchain interactions, α = 0.221(2) meV T −2/3 is a constant of proportionality related to (ℏ 2 /µ) 1/3 , and E0 = 0.92(2) meV is the energy required to create a pair of domain walls.The fit only includes data in the field polarized paramagnetic regime for 0.14 ≤ B ≤ 1.5 T; the zero field data were omitted from the fit as the magnetic structure is different from that above 0.14 T such that a different value of B0 would be needed.In both panels, the dashed line shows the energy of the lower boundary of the 2m1 continuum at l = 0.5 under the assumption that the quasiparticles do not interact.
is discussed in more detail in Sec.IV C. Fig. 7A also illustrates that the linear fit breaks down at low fields.This is expected, since at low field the number of spin flips ceases to be even approximately a good quantum number, and instead the weak confine-ment physics described in Sec.III A holds.In this low field regime, a fit was performed simultaneously to the band minima of all the modes, as shown in Fig. 7B, with the field dependence having a 2/3 power law form, as per (3), and the spacing between the energy levels being constrained by the zeros of the Airy function.The band minima were corrected for the effects of interchain dispersion by fitting the experimentally-extracted dispersion points to a parameterized three-dimensional (3D) dispersion relation of the form given in Appendix A 2 a, then setting the interchain hopping terms in those parameterizations to zero to correct for the small energy shifts due to interchain dispersion.The good agreement obtained between the field dependence of the bound state energies and the different expected power-law behaviours in the two field limits in Fig. 7 (solid lines) lends support to the proposal that, in CoNb 2 O 6 , the interchain couplings are sufficiently small relative to the dominant in-chain exchange, and the Ising character sufficiently strong, that both the weak and the strong spinon confinement regimes are realized experimentally via tuning of the applied field.The crossover region is approximately between 1.5 and 3.5 T. We note that, according to the Hamiltonian and parameters refined in Sec.IV, this crossover region is where the kinetic term for domain walls and the Zeeman term are comparable, such that the continuum approximation used in (2) for the weak confinement regime is no longer applicable.

IV. QUANTITATIVE DETERMINATION OF THE HAMILTONIAN
The Horace scans collected at high field have enabled us to further refine the microscopic model for this system beyond the minimal in-chain Hamiltonian proposed in [32].In this section, we will describe the refined inchain terms, as well as the fitting procedure, with details of the interchain interactions referred to Appendix A. There, we revise and extend a minimal model of the interchain interactions previously proposed to explain the dispersions in large transverse field [21], such that we can consistently reproduce the high field dispersions for field along both a and b directions.
The quantitative parameterization of the interchain dispersions obtained in Appendix A allows us to obtain effective 1D dispersion relations to which a singlechain Hamiltonian can be fit, which is done in Secs.IV A and IV B: Sec.IV A introduces the proposed single-chain Hamiltonian and Sec.IV B describes the fitting method.In Sec.IV C, we then demonstrate that this proposed Hamiltonian also quantitatively captures the behaviour of the magnetic excitations in CoNb 2 O 6 away from the fields at which the fit was performed.We also show that the proposed Hamiltonian can account for the spectrum observed in THz spectroscopy in [33].B) The Ising z-axis is at an angle γ to the c-axis in the ac plane.The absolute signs of λyz and angle γ cannot be determined from neutron scattering measurements, only their magnitude, so only one of the four possible orientations of the exchange ellipsoids compatible with the experiments is shown.The other three options are obtained by mirroring the exchange ellipsoid for a reference bond in a plane passing through its centre, parallel to the ab, bc, or ac crystallographic planes.The resulting orientation is then propagated via crystal symmetry operations onto all the other bonds, i.e. via the c-glide to obtain the other bonds on the same chain and via the b (or n glide) to relate bonds on the chain passing through the origin to bonds on the chain passing through the body-centre of the orthorhombic structural cell, as described in Fig. 14.

A. Proposed single-chain Hamiltonian
The proposed single-chain Hamiltonian is an extension of one recently proposed on symmetry grounds [32], and it is convenient to write it as where where j runs over all sites on a single chain.Here, xyz form an orthogonal right-handed coordinate system with y along the crystallographic b-axis and z defined to be the equilibrium spin direction in zero applied field, which for the Hamiltonian (9) coincides with the direction along which the spin components have the largest exchange, i.e. the Ising axis.We will show that H 1 is the minimal Hamiltonian required to qualitatively reproduce all key features of the the INS spectra seen in zero field, high field along a, and high purely-transverse field, while H 2 is needed to quantitatively capture all details of the spectra.
The first term in H 1 is the dominant Ising exchange.The second (λ S ) is a symmetric XY exchange term, which allows single spin flips to hop, and the third term (λ yz ) is an off-diagonal staggered exchange term which allows domain walls to hop.These three terms (together with H MF , defined below, for zero field and transverse fields below the critical point) can qualitatively account for all of the features seen in the INS spectrum, not only in field parallel to the a-axis presented here, across the full range of field values considered, but also in purely transverse field (∥ b) below [34] and above the critical field [32].In particular, the XY exchange is needed to account for the kinetic bound state seen in zero field (near l = −1 in Fig. 9A) [4], and for the large m 1 bandwidth seen in field parallel to the a-direction (Fig. 9E), as this term allows single spin flips to hop.The staggered off-diagonal exchange term λ yz is needed to account for the fact that domain walls can hop in zero applied field (seen in the fact that the bands around l = 0 in Fig. 9A are dispersive) [32].The key features of the nearest-neighbour exchange are graphically illustrated in Fig. 8: the major axis of the ellipsoids at the middle of each bond corresponds to the Ising exchange, the diameter in the perpendicular plane to the Ising direction illustrate the XY exchange λ S , and the staggered orientation of the ellipsoids is due to the staggered exchange λ yz .Finally, the last term in (10) is the Zeeman interaction with an applied magnetic field, where we have assumed that the g-tensor is diagonal in the xyz axes.
In order to quantitatively account for the full wavevector dependence of the data, however, other, subleading terms are needed, included under H 2 : the first term, λ AF , is a next-nearest-neighbour antiferromagnetic Ising term [4,32]  INS data in (A) zero field, (E) 8 T ∥ a and (I) 9 T ∥ b compared to dynamical correlations computed via ED for different models for a chain of 16 sites.Panel A is adapted from [4].White dots in B-D show experimentally-extracted dispersion points from the data in A, solid white lines in F-H and J-L show the analytic parameterization of the experimental dispersion relations with interchain terms set to 0. The models used were: the model 1D Hamiltonian in (9) with parameters in Table I (B, F, J); the model refined in [32] with gx and gz taken from Table I (C, G, K); the model used in [33] with gx and gz taken from Table I (D, H, L).In E-H, data below 5 meV (dashed line) are integrated over |h| < 1.5 and −1.1 < k < −0.9.Note that the shadow mode is the more intense mode here, as the range of transverse integration has been chosen in order to show the whole dispersion along l while minimizing artificial broadening due to interchain dispersion.Data above 5 meV are integrated over |h| < 1.5 and |k| < 1.5 and intensities in this region have been multiplied by 2 to make the m2 mode more clearly visible.In I-L, data below 4.1 meV (dashed line) are integrated over 0.8 < h < 1.2 and −0.1 < k < 0.1.Data above the dashed line have been integrated over −1 < h < 4 and |k| < 1 and intensities in this region have been multiplied by 4 to make the faint diffuse feature around l = −1 at about 5 meV more clearly visible.The ED calculations have been convolved with Gaussians of FWHM 0.066 meV, 0.16 meV and 0.3 meV in B-D, F-H and J-L respectively.In B-D, the intensity shown is S xx (l, ω), as the data in A have been corrected for the single ion magnetic form factor and the polarization factor under the assumption that S xx = S yy and S zz = 0 for inelastic scattering.In F-H and J-L, all components of the dynamical structure factor are included and the integration range, polarization factor and magnetic form factor have all been accounted for in the calculation.
bound state.The second term, λ xy AF , is an XY part to the next-nearest-neighbour term which accounts for the flattening of the bottom of the m 1 dispersion in high field along a [Fig.9E].The third term (λ A ) is an asymmetry between the XX and YY exchanges, which is needed to account for the position and bandwidth of the m 2 dispersion in large field along a; without this term, the m 2 mode is too low in energy.These last two terms are only of order 2-3% of the Ising term and the necessity for their inclusion in the parameterization can be seen by comparing Fig. 9E with G.The discrepancies illustrated in the latter between the empirical and calculated dispersions motivates our further refinement of the Hamiltonian with the inclusion of the λ A and λ xy AF terms which were fixed  to zero in [32].
For completeness, we note that, as mentioned in [32], two other nearest-neighbour exchange terms are symmetry allowed, λ xz J j S x j S z j+1 + S z j S x j+1 and λ xy J j (−1) j S x j S y j+1 + S y j S x j+1 .However, the definition of the axes used so far -that z is the direction of the equilibrium spin in zero field -places a constraint between these two terms, i.e., only one can vary independently.This is because each of these two terms on their own when added to (9) leads to a rotation of the zero-field equilibrium spin direction away from the z-axis, but for a given λ xz one can choose a corresponding λ xy of appropriate magnitude and sign such that, when both terms are present, the zero-field equilibrium spin direction is still along z.However, we find that allowing finite λ xz and λ xy with the above constraint does not measurably improve the agreement with the present experimental data, so in the following we assume λ xz = λ xy = 0.
Finally, H MF captures the effects of interchain couplings in a mean-field approximation, where in zero and low transverse field In fields above 0.14 T applied along the a-direction, the relevant form is instead given in Sec.IV C, and, at high field, this simplified form is no longer sufficient as excitations acquire a finite interchain dispersion at first order in the interchain couplings.Therefore, in this high field regime, we use the full explicit form for the relevant interchain exchanges proposed in (A4).
The refined parameter values are shown in Table I.The Hamiltonian shown above quantitatively reproduces the spectrum seen in zero field, high field ∥ a and high field ∥ b.Furthermore, it also reproduces data at low fields along b, as shown in [34] and at intermediate fields along a, and also accounts for previously published THz spectroscopy data in low transverse field, as we discuss later in Sec.IV C.

B. Fitting procedure
The fitting procedure used a global simultaneous fit to the dispersion relations corresponding to data taken in zero field, 8 T ∥ a and 9 T ∥ b, with the aim of arriving at a consistent description of all these different regimes within the same Hamiltonian.
First, cuts were taken through the data as a function of energy transfer at constant wavevector transfer.Dispersion points were obtained by fitting Gaussian peak shapes to these cuts.Many dispersion points were extracted from each data set (over 500 for the 8 T ∥ a data).Empirical dispersion relations were then fit to these dispersion points.For the high field data, these were 3D dispersion relations, as per (A2) and (A8) for field along a and b respectively.In zero field, the interchain hopping effects are negligible, due to the multi-spin nature of the bound states as well as the antiferromagnetic order pattern between chains, which suppresses interchain hoppings for the kinetic bound state (see Appendix B) so a 1D form was used.The single-chain Hamiltonian was then fit to these empirical dispersion relations with the interchain parameters set to zero.It was not possible to use the parameters derived from linear spin wave theory fits directly, since the predominantly 1D nature of the magnetic interactions together with the small effective spin S = 1/2 leads to strong quantum fluctuations that renormalize the dispersions even at the high magnetic fields investigated.In addition, we wanted to capture the various bound states, i.e. the zero field confinement bound states, as well as the high field ∥ a m 2 bound state; this is not possible in linear spin wave theory.Instead, exact diagonalization (ED) calculations on finite chains were used, with periodic boundary conditions at the ends of the chains.The fits used calculations on 12 sites as the best compromise between minimal finite size effects and the computation being quick enough to be carried out many times for fitting (for fitting to the kinetic bound state, 10 sites were used as this required more eigenstates to be found and is not strongly affected by finite size effects).Lanczos algorithms were used to diagonalize only the low energy subspace and speed up the calculation.There were 10 parameters in the fit but 13 pieces of information across these different data sets, so the fit was not underconstrained.In particular, the pieces of data used for the fit were the dispersions of the first two confinement bound states and the kinetic bound state from the zero field data [4], the m 1 and m 2 dispersions from the 8 T ∥ a data, and the magnon dispersion from the 9 T ∥ b data.For each of these data sets, the squared difference χ 2 between the ED at each momentum point and the empirical fitted dispersion relation was calculated, normalized by the uncertainty on the fitted dispersion relation as calculated from the covariance matrix of the fitted dispersion parameters, and summed; the fit minimized this total χ 2 .When calculating χ 2 for the high transverse field data, only the portion of the dispersion for l ≤ 0.5 was used, as at higher values of l, quasiparticle breakdown occurs [32,35], meaning that the dispersion relation ceases to be well defined.The optimization used a quasi-Newton algorithm and the initial parameters were varied in order to make sure that a global minimum was found.The very good one-to-one agreement obtained in the simultaneous fits of the ED to the empirical dispersion relations is shown in Fig. 10.
Uncertainties on the fitted parameters were estimated by varying the parameters of the empirical dispersion relations according to their covariance matrices.There are significant correlations between some of the parameters, especially between J and λ S and between λ S and g x .However, while a number of slightly different parameter sets give similar agreement with the features to which the single-chain Hamiltonian was fit, the final parameter set presented in Table I gives the best agreement not only with the features to which the Hamiltonian was fit, but also to other features in the spectra.These include the relative intensities of different features, the full bandwidth of the magnon in the high transverse field data, the energy and intensity of the faint diffuse feature around l = −1 at about 5 meV in the high transverse field data, the field dependence of the bandwidths of the m 1 and m 2 modes in field along a, and the spectra in low transverse field (see [34]).
A comparison between the data to which the Hamiltonian was fit and the spectrum calculated by exact diagonalization using this refined Hamiltonian is shown in the first two columns of Fig. 9.In these plots, the color indicates the measured/calculated scattering intensity, while the overlaid white dots/curves show the data points/dispersion relations that were being fit to.It can be seen that very good quantitative agreement is achieved.In these calculations, the plotted intensities are where α, β both run over x, y, z, f (Q) is the magnetic form factor, and Qα is the component along direction α of the unit vector parallel to the wavevector transfer Q.
The partial dynamical structure factors S αβ are where the sum is over all excited states |λ f ⟩ with energy E λ f relative to the ground state |GS⟩ and where the Fourier transformed spin is S α (Q) = r S α r e iQ•r , summing over all sites r.
The right hand columns of Fig. 9 contain comparisons to previous models that cannot quantitatively or qualitatively account for key features in the experimental data.These are discussed in Appendix C.This comparison provides evidence that all of the parameters are indeed needed in order to quantitatively capture all features in the data.

C. Comparison of the Hamiltonian to experiment at other fields
The refined Hamiltonian can also be compared to the results of INS data taken at fields other than those to which the fit was performed.The data agree significantly better with the calculation when using the fit with all terms included than when omitting any of the terms in H 2 .
Fig. 11 shows comparisons between ED calculations and the data at fields away from those used for the fits.Excellent quantitative agreement is found.In the calculations shown here, the component of interchain interactions parallel to the magnetization direction was taken into account in a mean field picture.For the field-polarized phase in field above 0.14 T along the adirection, the form of the interchain mean field in (12) does not hold because the ordering pattern of the chains changes, and instead the relevant form of the mean field term is where the interchain exchanges are defined in Appendix A 2 b.This was the form used in Fig. 11.In addition, in Fig. 12, we compare recently reported THz spectroscopy data (white dots) with predictions based on the Hamiltonian proposed here.The comparison is shown for fields up to 3.5 T, which we consider to be in the region where ED is sufficiently reliable, as the gap is still large (i.e., finite size effects are small), and interchain effects are also small.All the trends in the THz spectroscopy data are reproduced.The good agreement in the intensities can be seen visually by comparing Fig. 12 to Fig. 2c of [33] (not shown).There is good agreement over the whole field range probed.The dependence of the band maxima on field was found to discriminate more strongly between different models than the field dependence of the band minima.Note that the higher bound states cease to be well defined in ED once they overlap with the continuum; for this reason solid lines stop when they intersect the dashed line indicating the lower boundary of the 2m1 continuum in A. In B, the discrepancy between ED and experiment where the m1 and m2 ED curves touch is due to the fact that experimental values were extracted by fitting to the hopping form (A2) and points were not extracted in regions where different modes overlapped, as they do in Fig. 4C and D. The interchange in the gradient of the maxima of m1 and m2 modes with respect to field where they intersect is due to the exchange in character of the top of the mode: between approximately 2 T and 4 T, it is the top of the m2 mode that has single-spin-flip character around l = 1.

FIG. 12.
Dynamical structure factor at zero momentum transfer as a function of transverse field.Color indicates S xx (0, ω) as per ( 14), calculated using ED for a chain of 16 sites for the Hamiltonian in ( 9) and convolved with a Gaussian of FWHM 0.067 meV.White dots are data points extracted from the estimated local intensity maxima in the energy-dependent THz spectroscopy data presented in [33] at 1.5 K (Fig. 2c of that work).

V. CONCLUSIONS
In summary, we tuned the confinement potential between spinons in the ferromagnetic Ising chain material CoNb 2 O 6 by applying an external magnetic field along the crystallographic a direction, such that there was a large longitudinal (along Ising axis) field component.In the low field, weak confinement regime, we found a hierarchy of bound states, with their energy varying as field to the 2/3 power and the spacing between modes determined by the zeros of the Airy function, as expected in a picture of domain walls in a linear confining potential proportional to the strength of the longitudinal field.Upon increasing field, higher-order bound states increase more quickly in energy and progressively disappear from the spectrum such that in the limit of high field, in the strong confinement regime, we found only two bound states whose energies depend linearly on field.The higher energy of these two bound states is a dispersive two-spin-flip bound state which is stabilized by the proximity to the Ising limit.By performing a global fit to the full wavevector dependent spectrum observed in various fields, we also proposed a microscopic Hamiltonian including both in-chain and inter-chain interactions, down to 2% of the dominant Ising exchange.This Hamiltonian quantitatively reproduces the INS data obtained across a wide range of field conditions including zero field, high near-longitudinal field and high transverse field, as well as intermediate fields to which the Hamiltonian was not fit, leading to a fully-consistent description of the spin dynamics using a single set of exchange parameters.ACKNOWLEDGMENTS L.W. acknowledges very useful discussions with Michele Fava.We thank Alexander Chernyshev and César Gallegos for helpful comments and discussions.We acknowledge doctoral studentship funding from Lincoln College and the University of Oxford (L.W.), the Engineering and Physical Sciences Research Council (D.M.), and the University of Oxford Clarendon Scholarship Fund and NSERC of Canada (J.D.T).We acknowledge support from the Engineering and Physical Sciences Research Council grant numbers EP/H014934/1 (I.M.C. and R.C.) and GR/M47249/01 (D.P.).R.C. acknowledges support from the European Research Council under the European Union's Horizon 2020 research and innovation programme Grant Agreement Number 788814 (EQFT).R.C. also acknowledges support from the National Science Foundation under Grants No. NSF PHY-1748958 and PHY-2309135, and hospitality from KITP where part of this work was completed.The neutron scattering experiments at the ISIS Facility were supported by beamtime allocations from the Science and Technology Facilities Council.Figures 2 and 8 were made using the crystal visualization software vesta [36].Access to the data will be made available from Ref. [37].

Appendix A: Interchain dispersions in high field
In this Appendix, we present the characterization of the experimentally observed interchain dispersion relations at high field and their quantitative parameterization.

Characterization of the interchain dispersions in high field
INS data probing the dispersions in the interchain directions are plotted in Fig. 13.Panel A shows that in large field along a, no systematic dispersion bandwidth is detected along h within the experimental resolution used.Therefore, all data taken with field along a presented in the rest of this work were integrated across the entire range of h in the data, that is, over −1.5 < h < 1.5.Along k, however, there is visible dispersion and further, the bandwidth and even sign of the dispersion depends on the mode energy, i.e. on l.This is illustrated in Fig. 13B, where the lower energy and higher energy modes have k dispersions that are of opposite signs (white solid lines are fits as described in Appendix A 2 a).Note that the bottom mode is the shadow version of the primary m 1 mode wavevector-translated from l = 0 whereas the top mode is the primary m 1 mode at l = 1.This coupled kl dispersion is accounted for by the exchange pathways shown schematically in Fig. 13I and is discussed mathematically in the following subsection.
Panels D, E, G and H of Fig. 13, meanwhile, show the dispersions in the directions perpendicular to the chains in a transverse field, 9 T ∥ b.It can be seen that in this regime, there is interchain dispersion along both h and k, of similar magnitude to that seen along k for B ∥ a.However, while the dispersion along k remains of similar magnitude but switches sign when moving from l = 0 (Fig. 13E) to l = −1 (Fig. 13H), i.e. from the bottom to the top of the dispersion, similarly to what is seen for B ∥ a, the dispersion along h is strongly suppressed at l = −1 (Fig. 13G) compared to at l = 0 (Fig. 13D).The model proposed in [21] captures some, but not all, of these interchain dispersion effects.In particular it predicts that the interchain bandwidth is suppressed at higher energy (i.e.l = −1 compared to l = 0).This is indeed seen along h, but not along k.The same model also predicts that there will be very little interchain dispersion for near-longitudinal field, which is again seen along h but not along k; in fact, the dispersion bandwidth along k at l = 1 in near longitudinal field is approximately the same as in transverse field (the l = 0 transverse field k dispersion bandwidth is enhanced by the bonds in the ab-plane).We therefore expand and revise the interchain interaction model in [21] with the proposal that the bonds with a component along the adirection (the dash-dot light blue J 2 paths in Fig. 13J) are approximately Ising, whereas those in the bc plane [the solid and dashed light blue paths in Fig. 13I, i.e., J 1 (t 1 ) and J ′ 1 (t ′ 1 )] are approximately Heisenberg.

Parameterization of the dispersions relations in high field
Parameterization of the interchain dispersion relations was done by assuming that the dispersion of single spin flips in directions perpendicular to the chains could be quantitatively captured using a linear spin wave formalism, but with in-chain parameters which may be renormalized compared to those deduced from the exchange Hamiltonian.The linear spin wave formalism is expected to be asymptotically exact in the limit of high field (when the gap is much larger than the bandwidth).In this limit, the linear spin wave dispersion relation can be perturbatively expanded to become equivalent to a spin-flip quasiparticle hopping (or tight-binding) formalism.As the hopping formalism depends only on the exchange pathways, this same formalism can also be applied to quasiparticles composed of multiple spin-flips.In the case of field along the a-direction, we seek to parameterize the dispersion of bound states of various different numbers of spin flips within a single formalism, and therefore use a hopping formalism.For high transverse field, we seek to parameterize only a single magnon mode, but in a case where the gap is fairly small compared to the bandwidth such that a perturbative expansion in terms of a hopping To account for the coupled kl dispersion observed in high field ∥ a (see Fig. 13B), we propose a hopping model, shown schematically in Fig. 13I, where solid and dashed bonds indicate paths across which excitations could hop.This approach was used in order to be able to describe both the m 1 and the m 2 modes using the same formalism.The effects of the Ising parts of the interchain coupling are taken into account at a mean-field level (see Sec. IV C) which affects both the m 1 and m 2 states.This more phenomenological formalism is also able to capture the full 3D dispersions of higher bound states in lower fields.In Fig. 13I, the dark blue lines represent in-chain bonds (parameterized by t 3 and t 4 ), whose Hamiltonian is described in detail in (9), while the light blue lines represent inter-chain bonds (parameterized by t 1 and t ′ 1 ).The next-nearest-neighbour in-chain hopping (t 4 ) accounts for the flattening of the bottom of the dispersion [see Fig. 5A] while the diagonal interchain bonds (t ′ 1 ) explain why the dispersion bandwidth along k depends on l.
We can write this hopping model for the paths illustrated in Fig. 13I as a two sublattice model (a single zigzag chain with two sites per effective unit cell) with sublattices labeled α and β.The different Ising directions on the two chains in the structural unit cell can be neglected as we do not include interactions along the a direction since no dispersion along h is detected within the resolution of the experiment.This is consistent with the picture of the interchain couplings presented in the following subsection where the only coupling with a component along a, J 2 , is Ising like.This means that the hopping induced by J 2 is suppressed by a factor of sin 2 θ where θ is the angle between the local Ising axis and the local magnetization direction, and is small even for large fields along a.Therefore, in high field along a, we can approximate the Hamiltonian as being decoupled into bc planes with hopping interactions within planes and we write it as where c † R,α (c R,α ) creates (annihilates) a quasiparticle on the α sublattice in the unit cell with origin at R. The sum runs over all single-zigzag-chain unit cells, with primitive lattice vectors (a − b)/2, b and c (the projection of the primitive unit cell in the ab plane is shown by the shaded area in Fig. 14B).The dispersion relations in this two and the band minimum, E min = ω 0 + 2t 3 + 2t 4 (note that t 3 is negative for all bands as the dominant nearestneighbour interaction is ferromagnetic).It was found that, as the field decreases, the interchain dispersion decreases, which is consistent with the picture that upon lowering field, quasiparticles acquire a more pronounced multi-spin-flip character and the hopping of multi-spinflip excitations is suppressed as it is of higher order in the interchain couplings.The kinetic bound state remains dispersive down to 0.14 T, consistent with this being a single-spin-flip mode, but has unresolvable dispersion in zero field, consistent with the zero field magnetic structure, in which the interchain hopping of the single spin flip is suppressed, as shown in Appendix B. These trends were extracted from fits to the single orientation data, where the wavevector components k and l are coupled, and were confirmed by investigation of the Horace scan data set at 1.5 T.
We also note that the generic intensity from (A3) captures the general intensity dependence along the interchain directions as observed, e.g., in Fig. 13B.Together with the good agreement with the dispersions, this justifies a posteriori the use of a hopping model to parameterize the excitations in this regime.
b. Linear spin-wave theory parameterization of the 3D dispersions at high transverse field (∥ b) We here discuss the formalism used to describe the interchain dispersion in high transverse field (∥ b).
In the actual crystal structure of CoNb 2 O 6 , there are two zigzag chains per structural unit cell, shown in Fig. 14A (gray shading), one at the corner of the ab cell (filled circle) with Ising axis tilted at an angle +γ away from c towards a, and one chain in the center (open circle) with Ising axis tilted at an angle −γ.In order to make progress analytically, we work in a reference frame where the Ising axes of the central chains are rotated to match those of the chains in the corners, such that the unit cell halves (as shown in Fig. 14B), reducing the problem from a four sublattice to a two sublattice problem.This is possible because the interchain couplings are assumed to have a simplified form, being Heisenberg or Ising-like.In this frame, the Hamiltonian is where H 1,2 contain in-chain interactions defined in the (10,11).Here, R runs over all single zigzag chain unit cells, as in (A1).
We now solve this Hamiltonian in linear spin wave theory.Assuming spins polarized along the b-axis, after a Holstein-Primakoff transformation and a Fourier transform, the Hamiltonian in (A4) takes the form (up to quadratic order in magnon operators and omitting the constant term) where From this, we get the dynamical structure factor including the neutron polarization factors as where Q x,z are the components of Q projected along x and z respectively.In this expression, there are no mixed polarization terms as S xz = −S zx , and there are no cross-terms involving Fourier transformed spins at different wavevectors since these are zero by conservation of momentum.
To obtain intensities to compare with experiment, the delta functions in (A9) are replaced by Gaussians in energy to reflect the finite instrumental energy resolution and the intensity is multiplied by the squared spherical magnetic form factor for Co 2+ .The analytic calculations presented here were cross-checked against numerical calculations performed using SpinW [39].
The fact that the intensity contains terms with shifted wavevector Q + a * is a mathematical expression of the real space unit cell doubling (due to alternation of Ising axes) leading to Brillouin zone folding and a consequent shadow mode.This shadow mode is responsible for the additional weak scattering intensity at slightly higher energy than the main mode in Figs.13D and E. We see from (A10) that the intensity of the shadow mode relative to the primary mode is proportional to tan 2 γ.Note that this shadow mode is distinct from the shadow mode due to the zigzag of the chains, such that the full model in the fixed frame has a total of four dispersion relations, ω ± (Q) and ω ± (Q + a * ).
In order to extract a parameterization of the 3D dispersion in high transverse field, fits were done to the linear spin wave dispersion relations (A8).The fitted values for the interchain parameters are shown in Table II.We note that these parameters yield minima in the magnon dispersion relations at (1, ±q, 0) with q ∼ 0.39 which is compatible with the minima observed at 7 T ∥ b in Fig. 5 of [21], and can thus be related to the propagation vector q = 0.37 of the spontaneous incommensurate spin density wave order that sets in at the magnetic ordering temperature 2.95 K in zero field [11].To draw a direct comparison with interchain exchange parameters proposed in Ref. [21], the relevant value to compare to J 1 in [21] is J 1 + J ′ 1 in the present work.
Appendix B: Suppression of the interchain hopping of the kinetic bound state in zero field A surprising result found experimentally is that the kinetic bound state, the sharp mode near l = −1 at the top of the spectrum in zero field (Fig. 9A), has no measurable dispersion in either interchain direction, even though it is a single-spin-flip state, and generically one expects interchain dispersion at first order in the interchain couplings.We explain in this Appendix that this is due to the combination of the particular antiferromagnetic order pattern between chains and the form of the interchain couplings by introducing a simplified toy model.
As the interchain J 2 bonds in Fig. 13J are Ising like, and in zero field the spins are aligned along the Ising axis, the J 2 terms create no hopping along h.The J 2 bonds also do not contribute to a mean field, as the triangular lattice and the antiferromagnetic pattern mean that the effects of the two bonds cancel out.We therefore only need to consider the interchain bonds in the bc plane.In the simplest approximation, we consider straight (ζ = 0) ferromagnetic Ising chains along the c-direction, with an antiferromagnetic Heisenberg interaction between chains in the b direction, i.e., where the sum is over all magnetic sites.In the absence of any applied field, the chains order in an antiferromagnetic pattern along the b-direction (both in this simplified model and in the real material).Therefore, to calculate the linear spin wave spectrum, the quantization axis must be rotated by 180 • about the y-direction on alternating chains along the b-direction such that the ordered spin is along the +z direction on all chains.In this rotating frame, the quadratic spin wave Hamiltonian is This gives a dispersion relation ω(Q) = (J + J 1 ) 2 − J 2 1 cos 2 2πk.Now, J 1 ≪ J (experimentally, J 1 /J ≲ 2%), so this can be Taylor expanded, such that the bandwidth along k is proportional to J 2 1 /J, i.e. appears at second order in J 1 .This is not experimentally resolvable, which explains why no interchain dispersion is observed for the kinetic bound state in zero field.In contrast, in field ∥ a above 0.14 T, the spin components along a are all parallel, with the consequence that the interchain dispersion along k appears at first order in J 1 .This is experimentally resolvable, and indeed is clearly detected experimentally, with the kinetic bound state at 1.5 T having a bandwidth along k of 0.114(2) meV (not shown).

FIG. 1
FIG.1.A) Domain walls (solid lines) can separate apart at no energy cost in a ferromagnetic Ising chain.B) In a finite longitudinal field hz, there is an energy cost linear in the separation, as if the two domain walls interacted via a string tension λ (curly line between solid dots), which stabilizes confinement bound states.C) In the weak confinement regime (hz ≪ J), many bound states exist, which are coherent superpositions of states with the two domain walls separated by many sites.Solid lines are the Airy wavefunctions for the first three bound states as per (4) for parameters λµ/ℏ 2 = 0.072, relevant for the longitudinal mean field due to interchain interactions in CoNb 2 O 6[4].The domain wall separation, x, is in units of the spacing along the chain c/2.

FIG. 3 .
FIG.3.INS spectrum as a function of field from 0 to 2 T ∥ a, in the weak to intermediate confinement regimes, with field increasing from left to right and top to bottom.Color is the raw neutron scattering intensity (i.e., no background has been subtracted) on an arbitrary scale, collected using a high-resolution configuration (Ei = 4.30 meV) for a fixed sample orientation, with wavevectors projected along the chain direction.

4 .
Data at 2 T are shown in both configurations, [compare Fig. 3H and Fig. 4A].At 2 T (panel A), modes m 1 to m 4 (faint horizontal line near 3.7 meV) are clearly resolved, but, upon further increasing field, modes m 3 and m 4 become progressively fainter and above 6 T, in the strong confinement regime, only m 1 and m 2 remain (panels F to H).In the following we first discuss the strong confinement regime, then discuss in detail the intermediate field regime.

4 FIG. 4 .
FIG.4.INS spectrum as a function of field from 2 to 9 T ∥ a, covering the intermediate to strong confinement regimes, with field increasing from left to right and top to bottom.Color is the raw neutron scattering intensity on an arbitrary scale (different from Fig.3), collected using a high-coverage coarser-resolution configuration (Ei = 9.33 meV) for the same fixed sample orientation as in Fig.3, with wavevectors projected along the chain direction.Panel A is at the same field as the higher-resolution data in Fig.3H.

FIG. 5 .
FIG. 5. A) INS spectrum at 8 T, deep in the strong confinement regime, obtained by integrating a multi-angle scan over the transverse wavevector range |h| ≤ 1.5 and |k| ≤ 1.5.Intensities above the dotted line have been scaled up by a factor of 2 to make them more clearly visible.The incident energy used is the same as in Fig.4.B) Calculated intensity S(Q, ω) as defined in (13) for the single-chain Hamiltonian in(9), where interchain interactions are included in a meanfield approximation as per(15), calculated using ED on a chain of n = 16 (n/2 zigzag unit cells) sites with periodic boundary conditions.ED gives the spectrum for a discrete set of wavevectors spaced by ∆l = 2/n; for visualization purposes, the spectrum is plotted as constant over the interval ∆l around each discrete wavevector.The ED results have been convolved with a Gaussian in energy of FWHM 0.16 meV, representative of the experimental resolution in this energy range and obtained by a fit to the observed peak lineshapes.The white curves are the best fit single-chain dispersion relations with solid and dashed curves showing primary and shadow modes (obtained via l → l + 1 translation) respectively.The hatched region represents the 2m1 continuum obtained by assuming that quasiparticles do not interact.Note that m2 is below the 2m1 continuum at all wavevectors.The m2 and 2m1 intensities have been multiplied by factors of 2 and 40, respectively.

FIG. 6 .
FIG.6.A,B) Energies (blue symbols, left axes) and integrated intensities (red symbols, right axes) deduced from multi-Gaussian fits of the m3 (A) and m4 (B) mode at l = 0.5, extracted from cuts from data in the same configuration as in Fig.4.Solid lines are guides to the eye and the dashed lines show the estimated 2m1 continuum lower boundary at the same l-value.C-F) Energy scans showing the field dependence of the various modes (integrated over 0.4 < l < 0.6), solid lines are fits to multi-Gaussian peaks.Intensities are raw neutron scattering intensities on an arbitrary scale (the same for all panels).Note that a clear signal at the expected m3,4,5 energy is present even after crossing the 2m1 continuum lower boundary (dashed vertical line in panels C-F), which has been calculated using fitted dispersion relations for the m1 mode and assuming that m1 particles do not interact with each other.

FIG. 8 .
FIG.8.Graphical representation of the nearest neighbour part of the exchange Hamiltonian(9) projected onto the bc (A) and ac (B) planes for a single chain.The blue spheres represent cobalt atoms while the red ellipsoids at the midpoint of each nearest-neighbour bond represent the exchange matrix on that bond.The principal axes of the ellipsoids correspond to the principal axes of the exchange matrix, with the lengths of the principal axes proportional to the absolute values of the relevant eigenvalues of the exchange matrix.A) Consecutive bonds along the chain are symmetry-related by the c-glide of the crystal structure (mirror in the ac plane passing through the middle of each zigzag bond followed by translation by c/2), which leads to the staggered orientations of the ellipsoids between consecutive bonds along the chain; this corresponds to the staggered exchange term λyz in(10).B) The Ising z-axis is at an angle γ to the c-axis in the ac plane.The absolute signs of λyz and angle γ cannot be determined from neutron scattering measurements, only their magnitude, so only one of the four possible orientations of the exchange ellipsoids compatible with the experiments is shown.The other three options are obtained by mirroring the exchange ellipsoid for a reference bond in a plane passing through its centre, parallel to the ab, bc, or ac crystallographic planes.The resulting orientation is then propagated via crystal symmetry operations onto all the other bonds, i.e. via the c-glide to obtain the other bonds on the same chain and via the b (or n glide) to relate bonds on the chain passing through the origin to bonds on the chain passing through the body-centre of the orthorhombic structural cell, as described in Fig.14.
FIG. 9.INS data in (A) zero field, (E) 8 T ∥ a and (I) 9 T ∥ b compared to dynamical correlations computed via ED for different models for a chain of 16 sites.Panel A is adapted from[4].White dots in B-D show experimentally-extracted dispersion points from the data in A, solid white lines in F-H and J-L show the analytic parameterization of the experimental dispersion relations with interchain terms set to 0. The models used were: the model 1D Hamiltonian in(9) with parameters in Table I (B, F, J); the model refined in[32] with gx and gz taken from TableI(C, G, K); the model used in[33] with gx and gz taken from TableI(D, H, L).In E-H, data below 5 meV (dashed line) are integrated over |h| < 1.5 and −1.1 < k < −0.9.Note that the shadow mode is the more intense mode here, as the range of transverse integration has been chosen in order to show the whole dispersion along l while minimizing artificial broadening due to interchain dispersion.Data above 5 meV are integrated over |h| < 1.5 and |k| < 1.5 and intensities in this region have been multiplied by 2 to make the m2 mode more clearly visible.In I-L, data below 4.1 meV (dashed line) are integrated over 0.8 < h < 1.2 and −0.1 < k < 0.1.Data above the dashed line have been integrated over −1 < h < 4 and |k| < 1 and intensities in this region have been multiplied by 4 to make the faint diffuse feature around l = −1 at about 5 meV more clearly visible.The ED calculations have been convolved with Gaussians of FWHM 0.066 meV, 0.16 meV and 0.3 meV in B-D, F-H and J-L respectively.In B-D, the intensity shown is S xx (l, ω), as the data in A have been corrected for the single ion magnetic form factor and the polarization factor under the assumption that S xx = S yy and S zz = 0 for inelastic scattering.In F-H and J-L, all components of the dynamical structure factor are included and the integration range, polarization factor and magnetic form factor have all been accounted for in the calculation.

FIG. 10 .
FIG.10.Energies of dispersion points as calculated from ED on 16 sites compared to their values as calculated from the empirical dispersion relations.The one-to-one agreement obtained is excellent.

FIG. 11 .
FIG. 11.Band minima (A) and maxima (B) showing experimentally extracted points (symbols) corrected for interchain dispersion.The results of ED calculations on 16 sites (solid lines) are overlaid in the corresponding color, with effects of interchain couplings included as a mean-field correction.There is good agreement over the whole field range probed.The dependence of the band maxima on field was found to discriminate more strongly between different models than the field dependence of the band minima.Note that the higher bound states cease to be well defined in ED once they overlap with the continuum; for this reason solid lines stop when they intersect the dashed line indicating the lower boundary of the 2m1 continuum in A. In B, the discrepancy between ED and experiment where the m1 and m2 ED curves touch is due to the fact that experimental values were extracted by fitting to the hopping form (A2) and points were not extracted in regions where different modes overlapped, as they do in Fig.4C and D. The interchange in the gradient of the maxima of m1 and m2 modes with respect to field where they intersect is due to the exchange in character of the top of the mode: between approximately 2 T and 4 T, it is the top of the m2 mode that has single-spin-flip character around l = 1.

FIG. 13 .
FIG. 13.Interchain dispersions at 8 T ∥ a (A, B) and 9 T ∥ b (D, E, G, H). White curves represent best fit dispersion relations.Color is the raw neutron scattering intensity on an arbitrary scale.The incident energies Ei used were 9.33 meV (A, B), 2.14 meV (D, E) and 10.17 meV (G, H), and the overall intensity scale is different for each experimental configuration.C, F) Comparison of observed and calculated energies with best fit parameters for the 8 T ∥ a hopping dispersion (A2) and 9 T ∥ b linear spin wave theory dispersion respectively, with ω1,2(Q) = ω±(Q) and ω3(Q) = ω−(Q + a * ), with ω± as per (A8).I, J) Exchange interaction paths in the bc and ab planes, respectively.Labels t1,3,4 and t ′ 1 are hopping parameters for the dispersion model in (A2), while J1 and J2 refer to the Hamiltonian (A4).Dark/light colored lines indicate in-chain/inter-chain bonds, respectively.The gray shaded area indicates the structural unit cell shifted from the conventional P bcn unit cell such that the origin is at one Co 2+ site.The site labels α and β are used in the hopping calculations for the Hamiltonian in (A1) and in the linear spin wave theory calculations for the full interchain Hamiltonian in (A4).Filled/open circles represent spins with the local Ising axis at an angle ±γ to the c direction respectively.The experimental data were integrated in the transverse wavevector directions as follows: −1.05 ≤ k ≤ −0.95, −0.05 ≤ l ≤ 0.05 (A); −1.5 ≤ h ≤ 1.5, 0.95 ≤ l ≤ 1.05 (B); −0.1 ≤ k ≤ 0.1, −0.05 ≤ l ≤ 0.05 (D); 0.8 ≤ h ≤ 1.2, −0.05 ≤ l ≤ 0.05 (E); −0.1 ≤ k ≤ 0.1, −1.05 ≤ l ≤ −0.95 (G); 0.8 ≤ h ≤ 1.2, −1.05 ≤ l ≤ −0.95 (H).
FIG.14.A) ab plane of the crystal structure where filled/open circles represent Co 2+ ions with the local Ising axis at an angle ±γ away from c towards a.This alternation of the local Ising axis leads to two zigzag chains per orthorhombic unit cell (rectangular shaded area).These two chains are symmetryrelated by both a b glide (mirror in the1  4 ỹz plane followed by translation by b/2) and an n-glide (mirror in the xỹ1  4 plane followed by translation by (a + b)/2) of the P bcn space group of the crystal structure.Here xỹz define the fixed spin axes frame parallel to the orthorhombic abc crystal axes.Bottom diagrams show definition of a mathematically convenient spin axes frame xyz that rotates between the two chains such that z is always along the local Ising axis.B) In this rotating spin axes frame the two chains become equivalent and the unit cell halves (shaded parallelogram).
a † Q and b † Q are the magnon creation operators for the two single-chain sublattices.The matrix D(Q) has the formD(Q) =    A B C D * B * A D C C * D * A B D C * B * A    ,(A6)where A, B, C and D are functions of Q and are given Thus the Fourier transformed magnetization components are (in units of µ B )M x(Q) =g x S x (Q) cos γ + g z S z (Q + a * ) sin γ M z (Q) =g z S z (Q) cos γ − g x S x (Q + a * ) sin γ.

TABLE III .
Values of χ 2 corresponding to fixing some of the Hamiltonian (9) parameters to zero