Excitations of quantum Ising chain CoNb2O6 in low transverse field: quantitative description of bound states stabilized by off-diagonal exchange and applied field

We present experimental and theoretical evidence of novel bound state formation in the low transverse field ordered phase of the quasi-one-dimensional Ising-like material CoNb$_2$O$_6$. High resolution single crystal inelastic neutron scattering measurements observe that small transverse fields lead to a breakup of the spectrum into three parts, each evolving very differently upon increasing field. This can be naturally understood starting from the excitations of the ordered phase of the transverse field Ising model, domain wall quasiparticles (solitons). Here, the transverse field and a staggered off-diagonal exchange create one-soliton hopping terms with opposite signs. We show that this leads to a rich spectrum and a special field, when the strengths of the off-diagonal exchange and transverse field match, at which solitons become localized; the highest field investigated is very close to this special regime. We solve this case analytically and find three two-soliton continua, along with three novel bound states. Perturbing away from this novel localized limit, we find very good qualitative agreement with the experimental data. We also present calculations using exact diagonalization of a recently refined Hamiltonian model for CoNb$_2$O$_6$ and using diagonalization of the two-soliton subspace, both of which provide a quantitative agreement with the observed spectrum. The theoretical models qualitatively and quantitatively capture a variety of non-trivial features in the observed spectrum, providing insight into the underlying physics of bound state formation.


I. INTRODUCTION
The transverse field Ising chain (TFIC) is an important model in condensed matter physics because it displays the key paradigms of both a continuous quantum phase transition from an ordered phase to a quantum paramagnetic phase as a function of field, as well as, in the ordered phase, of fractionalization of local spin flips into pairs of domain wall quasiparticles (solitons) [1].The pure TFIC model can be mapped to non-interacting fermions [2,3], which in the ordered phase represent these solitons.However, a variety of different additional subleading terms in the spin Hamiltonian, such as a longitudinal field [4] or an XY exchange [5] can stabilize two-soliton bound states.Here we explore a regime where novel bound states can be stabilized by the interplay of applied transverse field and off-diagonal exchange.
The material CoNb 2 O 6 has been seen as a realization of TFIC physics for over a decade [5][6][7][8][9].Among the key experimental observations is the qualitative change in the nature of quasiparticles from domain walls in the ordered phase to coherently propagating spin flips in the highfield paramagnetic phase.Moreover, a fine structure of bound states was observed just below the critical transverse field, consistent with predictions for a universal E8 spectrum expected in the presence of a perturbing longitudinal field, which in this case arises from mean-field effects of the three dimensional magnetic order [5,9].The * These authors contributed equally to this work crystal structure is orthorhombic (space group P bcn), with Co 2+ ions with effective spin- 1  2 arranged in zigzag chains running along the c-axis, with dominant nearestneighbour ferromagnetic Ising exchange (see Fig. 1A).At the lowest temperatures, small three-dimensional interactions between chains stabilize a ground state with ferromagnetic ordering along the zigzag chains and with an antiferromagnetic pattern between chains [10][11][12][13][14].While the dominant magnetic physics in CoNb 2 O 6 can be captured by a TFIC Hamiltonian, additional terms in the Hamiltonian beyond the dominant Ising exchange are needed to explain various features of the spectrum [5,15].In particular, a staggered off-diagonal exchange term was recently proposed on symmetry grounds and shown to reproduce well the zero-field spectrum using density matrix renormalization group numerics [15].
In this work, we present high resolution single crystal inelastic neutron scattering (INS) data as a function of low to intermediate transverse field in the ordered phase.This regime has also been explored by THz spectroscopy [16], which probe the zone-centre (Q = 0) excitations.The INS data reveal a rich evolution of the magnetic spectrum with increasing field: the spectrum splits into three parts with each part behaving very differently.The top two parts are sharp modes, with the top mode becoming progressively flatter and the middle one progressively more dispersive in field, while the lowest energy part is a continuum where the intensity moves from bottom to top upon increasing field.
We seek to understand this rich behaviour in terms of a recently refined Hamiltonian model for CoNb 2 O 6 , which proposed all relevant additional exchange terms down to II.EXPERIMENTAL DETAILS Inelastic neutron scattering measurements of the magnetic excitations were performed on a large single crystal (6.76 g) of CoNb 2 O 6 grown using a floating-zone technique [18] and already used in previous INS experiments [5].The magnetic field was applied along the crystallographic b-direction, which is transverse to the local Ising axis of all the spins.The measurements were performed using the indirect geometry time-of-flight spectrometer OSIRIS at the ISIS facility.OSIRIS was operated with PG(002) analyzers to measure the inelastic scattering of neutrons with a fixed final energy of E f = 1.82 meV as a function of energy transfer and wavevector transfer in the horizontal (h0l) scattering plane.Throughout this paper, we express the wavevector transfer in the inelastic neutron scattering experiments as Q = (2πh/a, 0, 2πl/c) where (h, 0, l) are expressed in reciprocal lattice units of the orthorhombic structural unit cell, with lattice parameters a = 14.1337Å, b = 5.7019 Å and c = 5.0382 Å at 2.5 K [13].The sample was attached to the cold finger of a dilution refrigerator inside a vertical 7.5 T cryomagnet and measurements were taken at a temperature of 0.1 K.The average counting time at each field was around 7 hours.
For each field, two sample orientations were measured (c-axis oriented in the scattering plane at 25 • and 60 • with respect to the incident beam direction).Throughout this paper, the data panels presented are a combination of data from these two orientations, with the wavevector projected along the chain direction l as the physics considered is one-dimensional.The two orientations were chosen such that the projected l values covered a large part of the Brillouin zone along the chain direction.The INS data at one of the measured fields (2.5 T, in Fig. 2Q) was briefly reported in [15].
The data shown have had an estimate of the nonmagnetic background subtracted off, and have then been divided by the squared isotropic Co 2+ magnetic form factor f 2 (Q) and by the neutron polarization factor.The latter was calculated under the assumption that all in-elastic scattering is in the polarizations perpendicular to the Ising (z) axes and that the dynamical structure factor satisfies S xx (Q, ω) = S yy (Q, ω), an approximation which is found to be valid to a large extent for the model Hamiltonian (3) in the low transverse field regime.Here, where the sum extends over all excited states |λ f ⟩ of energy E λ f relative to the ground state |GS⟩ and where S x (Q) = j exp(iQ•r j )S x j , with j running over all sites.Under the above assumption, the wavevector dependence of the neutron polarization factor is Dividing the raw inelastic neutron scattering intensities by P(Q)f 2 (Q) then gives S xx (Q, ω) up to an overall scale factor.Eq. ( 2) is appropriate for the experimentally observed zero-field magnetic structure of CoNb 2 O 6 and takes into account the two different chains per crystallographic unit cell with Ising directions at an angle of ±γ to the c-direction in the ac-plane.We have taken γ to be 30 • [13].

III. EVOLUTION OF THE MAGNETIC EXCITATIONS WITH APPLIED FIELD
In this section, we first introduce the model Hamiltonian and relate this to the zero field spectrum, introducing the concept of two-soliton states.In applied field, the spectrum splits into three components with different evolution in field.We show in the following sections that this rich behaviour can be naturally understood in a picture of solitons with dispersions tuned by the transverse field.

A. Model Hamiltonian
We use the single-chain Hamiltonian model for CoNb 2 O 6 recently refined in [17].It is convenient to write this in three parts: where and and where z is the Ising direction (in the crystallographic ac plane, defined as the direction of the moments in zero applied field), y is parallel to the crystallographic b-direction and x completes a right-handed coordinate system.This Hamiltonian, with parameter values in Table I, is a refinement of the minimal model proposed in Ref. [15] and was recently deduced using a simultaneous fit to the spectrum in zero field, large transverse field, and large near-longitudinal field [17].One can regard H 1 as the minimal Hamiltonian needed to qualitatively reproduce all key features of the excitation spectrum, while the terms in H 2 are sub-leading and are added in order to achieve quantitative agreement.This applies to both the data in Ref. [17] as well as data in low transverse field presented here.Finally, H 3 captures the effects of the weak interchain interactions at a mean-field level.The given form with a constant λ MF > 0 applies throughout the field range explored here (0 to 2.5 T ∥ b) as the magnetic order pattern between chains does not change in this field range [19] and as, for this magnetic structure, all the interchain interactions that have a net contribution to the mean field are Heisenberg-like [17].The different signs in front of the S y j and S z j terms reflect the fact that the S y -components of spins on neighbouring chains are parallel (polarized by the applied field h y ), whereas the S z components are (spontaneously) aligned antiparallel by the antiferromagnetic interchain interactions.
The dominant term in the model Hamiltonian is the first term in H 1 , the nearest neighbour ferromagnetic Ising exchange.The second term is a nearest neighbour ferromagnetic XY exchange term (λ S ), which causes single spin flips to hop.The third term is a staggered off-diagonal exchange term (λ yz ) which causes solitons to hop with a sign that alternates along the legs of the zigzag chain.The transverse field term (h y = g y µ B B y ) flips single spins; for a fixed number of domain walls, this is equivalent to soliton hopping by one site.We will show that the competition between these two one-soliton hopping terms leads to a rich field-dependent spectrum.
The spectrum in zero field (Fig. 2A data, 2B calculation) can be qualitatively understood in terms of the minimal Hamiltonian H 1 + H 3 .A spin-flip neutron scattering event creates a pair of solitons (domain walls) and FIG. 2. Evolution of the INS spectrum and calculated spectrum for the full Hamiltonian in (3), as a function of applied magnetic field, increasing from top to bottom.Columns left to right: Inelastic neutron scattering data (A adapted from [5]); S xx as calculated by ED on 16 sites (8 unit cells) with periodic boundary conditions; S xx as calculated by diagonalization of the two-soliton subspace (see Sec. IV B) on 100 sites (50 unit cells) with periodic boundary conditions; S xx as calculated analytically by perturbing about the localized limit for λMF=0 (see Sec. V).In the right-most column, the hatched patches represent the two-soliton continua; in this approximation, the continua have no scattering intensity except where the bound states overlap them and create a resonance.In each panel, color indicates S xx , as defined in (1) and further normalized by the total number of sites, on a linear scale indicated by the colorbar.The calculations have been convolved with a Gaussian of FWHM 0.067 meV to mimic the estimated experimental resolution effects and calculated intensities are shown in absolute units of meV −1 .The intensities for the data in panels E, I, M, Q, have been multiplied by a common scale factor to bring them visually into agreement with the corresponding calculations in column 2. The data in A come from a different experiment so a separate scale factor was used for those intensities to bring them into agreement with those in panel B.   4)-( 6), from Ref. [17] for the pure Ising chain, the energy is independent of the separation between these solitons [see Fig. 1C].In the presence of the staggered off-diagonal exchange (λ yz ), the solitons become dispersive [15], resulting in a continuum of scattering in energy-momentum space, covering a large energy extent near l = 0.A longitudinal interchain mean field (H 3 ) acts as an effective linear potential confining the solitons into a series of bound states, as seen in Fig. 2A near l = 0.The sharp mode in the data near l = 1 is a two-soliton kinetic bound state stabilized by the XY (λ S ) exchange [5].
The terms in H 2 do not change the qualitative content of the spectrum but are important for quantitative agreement with the experimental data.The first term in H 2 is an antisymmetric diagonal nearest neighbour exchange (λ A ).The second and third terms are a next nearest neighbour antiferromagnetic XXZ exchange (λ AF and λ xy AF respectively).The second term is needed to account for the energy of the kinetic bound state near l = 1 [5], while the first and third are needed to explain the details of the dispersions seen in very high field [17].Fig. 2B demonstrates that very good quantitative agreement has been achieved between the experimental zero-field data and exact diagonalization calculations using the full Hamiltonian (3).

B. Spectrum in small to intermediate transverse field
The key feature of the evolution of the spectrum as a function of field, shown in the first column of Fig. 2, is the break-up of the observable spectrum into three parts, each evolving very differently upon increasing field.The top part, which evolves out of the zero-field high-energy kinetic bound state, is a sharp mode that becomes progressively flatter upon increasing field and spreads out over the whole Brillouin zone.In contrast, the middle part is dominated by a sharp mode that becomes progressively more dispersive upon increasing field and appears to trade intensity with the top mode.The lowest energy part is dominated by a continuum spread, with intensity moving from bottom to top upon increasing field.All features and trends in the INS data are quantitatively captured by exact diagonalization (ED) calculations using the full Hamiltonian in (3) with the expectation values ⟨S y ⟩ and ⟨S z ⟩ in the mean-field term H 3 calculated selfconsistently; these calculations are shown in the second column of Fig. 2.
The breakup of the spectrum in field is clearly illustrated in Fig. 2I at 1 T. The lowest energy part of the spectrum centred around 1.4 meV shows a set of excitations extending over a broad energy range with clear sharp modes visible both at the bottom and the top of this range.Those excitations are clearly separated from another set of states centred around 2.25 meV with a clear sharp mode near 2.5 meV.At higher energies still, there is the vestige of the zero-field kinetic bound state near l = 1, now clearly separated from the rest of the spectrum.All the observed features, both dispersions and wavevector-dependence of intensities, are well captured by the ED calculations in Fig. 2J.
The 2.5 T data (Fig. 2Q) shows even more contrasting behaviour between the different parts of the spectrum.The sharp mode at the top of the low energy continuum now extends all the way from the Brillouin zone center to l ≈ 0.6 and has gained in intensity compared to the continuum below it.In the middle energy region the sharp mode has become strongly dispersive, with the middle continuum losing nearly all its scattering intensity, and the top sharp mode has become almost entirely flat and spread out over almost all the Brillouin zone.Again, all these features are well reproduced in Fig. 2R.The spectra at 0.5 T (Figs. 2E and F) and 1.5 T (Figs. 2M and N) interpolate between 0, 1 and 2.5 T and show the gradual evolution of the spectrum.
The ED calculations quantitatively capture every feature and trend described above.We stress that the parameters used in this calculation were not fit to the finite transverse field data presented here, but are fixed to the values proposed in [17].This excellent agreement between data and calculation gives further support to the Hamiltonian proposed in [17] and motivates our search for a physical picture of the excitations.In the following sections, we will introduce a picture of solitons on the zigzag chains and show that the breaking up of the spectrum and very different evolution of the different parts in field can be captured quantitatively and understood phenomenologically in terms of solitons hopping and bound state formation.

IV. TWO-SOLITON MODEL
In (4) to (6), all λ terms are ≪ 1 such that the dominant term is the ferromagnetic Ising term.This means that it is sufficient for our purposes to consider two-soliton excitations, and neglect mixing with fouror-more-soliton excitations, since those occur at much higher energy.More systematically, we may write the Hamiltonian as with H Ising denoting the Ising Hamiltonian, and V containing all other terms in H.We now treat V as a perturbation of order δ, and define a Schrieffer-Wolff transformation, i.e., that H ′ conserves the number of solitons (domain walls).By expanding S = S 1 δ + S 2 δ 2 + ..., and enforcing (7) up to order δ n with n ∈ Z + , we obtain a systematic perturbative series for the effective Hamiltonian within subspaces with a fixed number of solitons.Up to lowest order in δ, we arrive at with V conserv denoting the soliton number conserving terms in V, and P i standing for the projector to the sector with exactly i solitons.We have verified that the second order terms, ∼ δ 2 , are negligible compared to this leading contribution.
Relying on these insights, we now start by considering the effect of the effective Hamiltonian ( 8), first on a single soliton, and then within the two-soliton sector.We focus on the minimal Hamiltonian H 1 + H 3 , yielding a good qualitative understanding of the spectrum, and leave the discussion of H 2 to the Appendices.We further assume that the most important contribution from the mean field interchain coupling H 3 is a z magnetic field, For convenience, we define the unconventional raising and lowering operators, which obey the usual commutation relations.Note that this is equivalent to performing the calculations in a spin basis rotated by π/2 around the z axis, obtained via the canonical transformation S x j → −S y j and S y j → S x j in H.

A. Action of the Hamiltonian on a single soliton
In this section we consider the spectrum of deconfined solitons under the Hamiltonian H 1 projected to the single soliton sector.The confining mean field H 3 will be introduced in Sec.IV B, where we discuss the spectrum in the two-soliton sector.
Let us define the "left" soliton state, a single domain wall at the link (j − 1, j), separating up spins on the left and down spins on the right, Here, the arrows indicate the eigenstates of S z j with eigenvalues ±1/2.We now consider the action of the Hamiltonian H 1 on this state, term by term.
• Ising exchange: The single soliton state |j⟩ L is an eigenstate, with excitation energy ϵ 0 = J/2 above the ground state.
• XY exchange λ S : This term flips two adjacent spins in opposite directions.Acting on |j⟩ L , it creates two additional domain walls by flipping the spins at sites j − 1 and j, and can therefore be dropped.
• transverse field h y : This term flips a single spin, where we have used the unconventional raising and lowering operators (9).To conserve the number of solitons when acting on |j⟩ L , the flipped spin must be either at j or j − 1, Therefore, the transverse field gives rise to a nearest neighbour hopping term for the domain wall.
• staggered off-diagonal exchange λ yz : Similarly to the transverse field h y , this term results in a single spin flip, Here, the operator S + j + S − j can flip the spin at site j if and only if the spins on sites j + 1 and j − 1 are opposite, due to the factor S z j+1 − S z j−1 .Therefore, V yz yields spin flip processes confined to domain walls.We find a hopping term similar to the effect of the field h y , but with a different sign structure across links.
To obtain the spectrum of this hopping Hamiltonian, we write a Schrödinger equation for the soliton.We define the state where ψ L (j) = L ⟨j|ψ⟩ is the wavefunction of the soliton.Using the expressions derived for P 1 H 1 |j⟩ L above, the Schrödinger equation, This equation describes a staggered hopping of solitons, with hopping amplitudes alternating on even/odd bonds.Since H is invariant under translations by two lattice sites, we can use the following Bloch ansatz, where k = 2πl/c is the soliton momentum.In the following, we will interchangeably use the symbols k and l when referring to momentum along the chain direction, with the only difference that k is in absolute units whereas l is in reciprocal lattice units.In the above equation, the coefficients ψ Lσ differentiate between the even and odd sublattices, reflecting the two-site unit cell of the Hamiltonian.From now on, we will reserve the index p for labelling the unit cells, whereas j will be used as a label of lattice sites.With this convention, the Schrödinger equation reduces to yielding a pair of bands, ω ± , with bonding / anti-bonding orbitals Importantly, the dispersion vanishes if h + = 0 or h − = 0.In these limits, the hopping amplitude vanishes either on odd or on even bonds, and the domain wall can only move between two sites, giving rise to flat localized bands, as previously noted in [16].This localized limit will serve as a convenient starting point for perturbative considerations in Sec.V, allowing us to obtain a simple qualitative picture for the evolution of the INS spectrum with magnetic field h y .Above we considered a single "left" soliton, describing a domain wall with up spins on the left and down spins on the right.Another type of domain wall excitation is a "right" soliton, separating a domain of down spins on the left from up spins on the right, The arguments described above can be repeated for right solitons, with the only difference being that the hopping induced by λ yz is of opposite sign compared to the case of left solitons.This interchanges the hopping amplitudes h + and h − in the Schrödinger equations, but leaves the dispersion (12) unaltered.Therefore, both solitons become localized at the same critical magnetic field h y .

B. Solution of the Hamiltonian in the two-soliton subspace
We now turn to the spectrum of the effective Hamiltonian within the two-soliton subspace.In Sec.IV A we obtained two distinct soliton dispersions, ω ± , describing bonding/anti-bonding orbitals.Therefore, we expect three continua in the two-soliton subspace, arising from the pairings (ω + , ω + ), (ω + , ω − ) and (ω − , ω − ).However, the solitons interact, due both to hard-core repulsion and to nearest neighbour soliton-soliton interactions encoded in H 1 .Moreover, the full Hamiltonian also includes a confining z magnetic field, H 3 , yielding an attractive interaction between the two solitons in a pair.Below, we take into account all of these effects by considering H 1 + H 3 projected to the relevant subspace with P 2 and show that a full understanding of the spectrum cannot neglect these interactions.
Assuming h z > 0, the relevant low energy excitations correspond to a single domain of down spins inserted into a background of up spins.Therefore, it is convenient to use the following basis, with a "left" soliton on the left and a "right" soliton on the right, with j L < j R .Similarly to the procedure followed in Sec.IV A, we can derive a Schrödinger equation within the two-soliton subspace by considering the effect of P 2 H 1 and H 3 on the basis states.
The transverse field h y and staggered off-diagonal exchange λ yz again give rise to hopping terms for the left and right solitons, with a correction term arising for nearest neighbor solitons j R = j L + 1 due to hard core repulsion, Here ′ denotes a restricted summation constrained to valid basis states, by dropping the unphysical terms |j L + 1, j R ⟩ and |j L , j R − 1⟩ for neighboring solitons j R = j L + 1.
Besides these familiar terms, two new types of contributions arise compared to the single soliton case.The XY exchange term, gives rise to a nearest neighbor interaction term between solitons, This term shifts the center of mass coordinate of the soliton pair by one lattice site, without changing the relative coordinate j R −j L .Finally, the magnetic field h z leads to an attractive potential between the left and right soliton, Based on these relations, we construct the two-particle Schrödinger equation for the wavefunction Ψ(j L , j R ) defined through Relying on translational invariance for the center of mass coordinate, it is convenient to write The two-soliton Schrödinger equation derived above yields a spectrum consisting of three continua and three bound states across a wide range of parameters, see Appendix C. As mentioned above, the origin of the three continua can be understood as due to the three different ways of combining the two bands of solitons into twosoliton continua.The origin of the bound states, which we term ε bound states to avoid confusion, will be explored in the following section, Sec.V.
Before turning to the detailed study of the ε bound states, we conclude this section by deriving a formula for the INS spectrum within the two-soliton model, showing that the dominant contribution stems from the three bound states.To this end, we approximate the ground state |GS⟩ appearing in the dynamical structure factor S xx , Eq. ( 1), as the ground state of the Ising Hamiltonian, |GS⟩ ≈ |... ↑↑↑ ...⟩.Acting with the spin operator S x (k) creates a soliton pair.Using the unconventional raising and lowering operators S ± j = S y j ∓ iS x j , we can write the resulting state as where k = 2πl/c is the soliton pair center of mass momentum.By substituting the eigenstates |λ f ⟩ with the solutions of the two-soliton Schrödinger equation constructed above, we arrive at the overlaps Therefore, the ε bound states, which we will show to have a large weight on the configurations with nearest neighbor domain walls, give the dominant contribution to the dynamical spin structure factor (1).
The INS intensity as calculated above is shown in the third column of Fig. 2. The calculation uses the full Hamiltonian with the same parameters as in the ED calculations, except that the spin vector expectation value ⟨S⟩ = (0, 0, 1/2) is assumed fixed, rather than using a self-consistent value.This is a good approximation since even at 2.5 T, the self-consistent value as calculated by ED is ⟨S⟩ = (0, −0.150, 0.473).The agreement between the observed spectrum and the model is still quantitative -all features and trends are captured -but not quite as strong as for the exact diagonalization calculation.For instance, there is a small overall energy shift, most visible by comparing Figs.2R and S; in the latter, energies are shifted to lower values.However all key features are well reproduced at all measured fields.
To gain more insight into the structure and magnetic field dependence of this INS signal, we examine the ε bound states in the next section, by relying on a perturbative argument around the localized limit h − = 0.

V. THE LOCALIZED LIMIT
As derived in Sec.IV A, the staggered off-diagonal exchange λ yz and the transverse field h y in the leading order Hamiltonian H 1 lead to hopping terms of opposite sign for the solitons.A particularly interesting situation arises when these terms are matched, such that h − = 0, resulting in localized single solitons.This localized limit serves as a convenient starting point for perturbative considerations, shedding light on the structure and magnetic field dependence of the ε bound states, as well as the three two-soliton continua.First, in Sec.V A, we set h − = 0, and study the two-soliton spectrum, in particular, the nature of the two-soliton bound states.We then examine the effect of a small non-zero delocalizing term h − in Sec.V B. Predictions for the evolution of the INS spectra with decreasing transverse field h y , and comparisons of these to the experimental data, are discussed in Sec.V C. For most of this section, we focus on the leading order Hamiltonian H 1 , with a brief comment about the mean field H 3 at the end of the section....
3. Soliton hopping in the localized limit.A) In the case of well separated solitons, they can hop independently, with each soliton hopping between two sites with matrix element h+.B) When solitons are in adjacent unit cells, hopping is constrained by hard-core repulsion between solitons.

A. Localized limit h− = 0
For simplicity, we start by setting λ S = 0 as well, and only keep the staggered off-diagonal exchange λ yz and the transverse field h y .We will discuss the effect of the nearest neighbor exchange λ S later.Under these simplifications, a left soliton can hop between the two sites of a unit cell p, 2p ↔ 2p + 1, with rate h − , whereas it hops between neighboring unit cells p − 1 and p, through sites 2p − 1 ↔ 2p, with rate h + , see Fig. 3A.For h − = 0, we obtain the following eigenstates with energies ω ± , symmetric and antisymmetric under the inversion exchanging the even and odd sublattices, respectively.For a right soliton, the role of h − and h + is interchanged, leading to symmetric / antisymmetric eigenstates localized within a single unit cell p, Relying on these observations, we can construct the localized two-soliton eigenstates by considering a left soliton confined to sites 2p − 1 and 2p, and a right soliton on 2p ′ and 2p ′ + 1.For p ′ > p, the solitons do not interact, and we obtain the eigenstates The eigenvalue ω +− = J is doubly degenerate, and the eigenstates were chosen to be symmetric / antisymmetric under inversion.We now construct delocalized eigenstates with a well defined center of mass momentum k as follows with N denoting the number of unit cells, and n ≥ 1.These eigenstates correspond to the three two-soliton continua arising from the different pairing of bonding / anti-bonding orbitals.In the localized limit considered here, we obtain highly degenerate flat bands at energies J ± 2h + and J, reflected by the free index n standing for the relative coordinate between the left and right solitons.
Placing the left soliton to sites 2p − 1 and 2p, and the right soliton to 2p ′ and 2p ′ +1 with p = p ′ gives rise to interaction through hard core repulsion, see Fig. 3B.In this case, the eigenstates can be obtained by diagonalizing a 3 × 3 matrix acting on the three allowed configurations, yielding with energies ε ± = J ± √ 2h + and ε 0 = J.The corresponding momentum eigenstates form non-degenerate flat bands (Fig. 4A) given by These states are the origin of the ε bound states found through the numerical solution of the two-soliton Schrödinger equation in Sec.IV.We now consider the effect of a weak symmetric exchange λ S , while keeping h − = 0.As discussed in Sec.IV, this term only affects nearest neighbor solitons by shifting the center of mass coordinate.Therefore, the three continua, (15), with solitons residing in different unit cells are not affected.In contrast, the symmetric exchange acts non-trivially on the bound states ( 16), inducing an energy shift calculated perturbatively as The full spectrum of the localized limit h − = 0, with weak symmetric exchange λ S is illustrated in Fig. 4B, showing the three highly degenerate flat continua, and the three dispersive ε bound states, which become delocalized by λ S .Note that in CoNb 2 O 6 , λ S is large enough that the dispersion causes the ε 0 and ε + bands to cross, leading to band inversion.This effect is discussed quantitatively in Appendix D and illustrated in Fig. 4C.The band inversion further suppresses the dispersion of the top mode, as well as mixing the character of the two bands.

B. Effects of weak delocalizing hopping h−
As shown above, the localized limit h − = 0 provides remarkable insight into the structure of the two-soliton spectrum obtained near to this limit.We can gain a more detailed understanding of the evolution of the spectrum with decreasing transverse field by considering the effect of a weak delocalizing hopping term h − perturbatively.The first important effect of h − ̸ = 0 is lifting the high degeneracy of the three continua and broadening these bands.Secondly, a finite h − can mix the ε bound states constructed in the previous section with the continua in k regions where they overlap in energy.
We first explore the first effect, by focusing on the continua, and applying degenerate perturbation theory within each band separately.We note that for weak h − and λ S , the bottom and top bands remain well separated from the ε bound states.Assuming also λ S ≫ h − , the middle continuum overlaps with the middle bound state only in the vicinity of k = π/c, see Fig. 4D.Therefore, our treatment of neglecting the hybridization between the continua and the ε bound states is justified in the limit of weak couplings h + ≫ λ S ≫ h − , apart from the case of the middle band in the vicinity of k = π/c.While the experimental parameters lie outside of this well controlled region, we will demonstrate below that a first order perturbative expansion grants valuable insight into the evolution of the INS spectra for the whole range of applied transverse fields.
Denoting the hopping term by V h− , we find that the matrix elements between the single soliton eigenstates of the localized limit are Up to first order in perturbation theory, the energy shifts of the three continua can be evaluated by considering the matrix elements of V h− between the two-soliton eigenstates of the localized limit, (15), within each band separately.Relying on the relations (18), we obtain for the top and bottom bands For a fixed center of mass momentum k, this equation corresponds to an effective nearest neighbor hopping Hamiltonian for the relative coordinate n, with hopping amplitude ±h − cos(kc/2), subject to the hard core constraint n > 0. Therefore, for the bottom and top bands we find the spectrum, with approximate unnormalized eigenstates Here, k is the total momentum of the soliton pair and q is the relative momentum of the two solitons, and the factor sin(qnc/2) reflects hard core repulsion n > 0. Thus, the weak hopping h − broadens the bands, most strongly around k = 0, but the high degeneracy still persists at k = π/c, see Fig. 4D.Turning to the middle continuum, we have to calculate four types of matrix elements between the two-soliton states |n, k⟩ 0 ± .We obtain We can now calculate the broadening of this band by diagonalizing this matrix within a given total momentum sector k.We find that the eigenstates in the middle band remain at least twofold degenerate everywhere, yielding the spectrum with approximate unnormalized eigenstates Dispersions and intensities as a function of wavevector and energy at different stages of perturbation away from the purely localized limit.The x-axis is momentum in reciprocal lattice units, l = kc/(2π), vertical axis is energy relative to J.In each panel, color indicates S xx , as defined in (1) and further normalized by the total number of sites, on a linear scale indicated by the colorbar; black lines represent ε bound states while gray hatched patches represent continua.The calculated intensities have been convolved with a Gaussian of FWHM 0.033J and are shown in absolute units of 1/J.In all panels, h+ = 0.2J.
In A-D) the calculation uses the perturbative regime of Sec.V. A) Bound states and continua are both dispersionless for h− = 0 and λS = 0. B) For finite λS, bound states become dispersive but continua remain flat and highly degenerate.C) Band inversion occurs between the top two bound states ε+ and ε0 for JλS > 2 √ 2h+/3 (see Appendix D), resulting in modes labelled ε1,2, such that ε1 inherits the structure of the state ε0 close to l = 0 and l = 1 while retaining the character of ε+ elsewhere, with ε2 showing the opposite behavior.This is reflected by the change in intensity distribution from B to C near l = 0, 1. D) Continua become dispersive for finite h−.E) Same as D, but calculated using the full two-soliton model in Sec.IV, which predicts hybridization between the lowest two bound states, ε2 and ε−, and the continua they overlap.
with ξ(n) = ± for n even/odd, and q again standing for the relative momentum.Thus, the middle band remains highly degenerate around k = 0, but it is broadened away from this point, most strongly around k = π/c.
We note that these dispersions of the continua can be understood based on the single-soliton spectrum as given in (12).In the limit of small h − , (12) becomes The three continua constructed above correspond to the three different types of soliton pairs.Denoting the individual solitons' momenta by k 1 and k 2 , we obtain the energies with total wavevector k = k 1 + k 2 , and relative momentum q = k 1 − k 2 , in accordance with the expressions derived above.
The second important effect of the hopping h − is mixing the ε bound states with the continua where they overlap in energy.In these regions, the originally bound soli-ton pair can become delocalized, broadening out the INS signal, as shown in Fig. 4E.

C. Comparison with INS data
We conclude this section by describing the INS intensity predicted by these perturbative arguments.In the localized limit h − = 0, the continua (15) have no overlap with nearest neighbor soliton pairs |j, j + 1⟩ so do not contribute to the overlap (14) and to the resulting INS spectrum.The ε bound states (16), on the other hand, have the property that following from comparing ( 16) to (13).Therefore, the different bound states contribute differently to the INS spectrum S xx (Q, ω).This structure is inherited by the bound states away from the limit h − = 0.
In the absence of band inversion (Jλ S < 2 √ 2h + /3), the three bound states can be labeled by ε ±,0 , such that the bottom and top modes ε ± yield strong signals in the vicinity of l = 0 and are suppressed near l = 1 [note that l = kc/(2π)], whereas the middle band ε 0 behaves in the opposite way, see Fig. 4B.If λ S is large enough for band inversion to occur, the top and middle bands acquire new labels ε 1,2 , with ε 1 retaining the structure of ε + around l = 0.5 but inheriting the character of ε 0 around l = 0 and l = 1, and the reverse holding for ε 2 .This results in a transfer of intensity between the top two bound states in the vicinity of l = 0 and l = 1.That is, the top mode is strong at l = 1 and weak at l = 0, and this is opposite to the middle band (Fig. 4C).
For small but finite h − , the lowest ε bound state mode ε − , which is strong around l = 0 in the limit h − = 0, is pushed into the bottom continuum due to the strong broadening of the continuum with h − .The mixing between these states leads to a signal smeared out across a larger range of energies.Similarly, the middle bound state mode mixes with the middle continuum around l = 0.5, smearing and eventually almost completely washing out the signal, see Fig. 4E.However, for λ S large enough to lead to band inversion, the top mode does not hybridize with the continuum around l = 1, even for large h − .This is because the upper continuum consists of states that are even under exchange of the even and odd sublattices, whereas the bound state ε 0 is odd, shedding light on the remarkably sharp INS spectrum of the top state around l = 1, even far from the localized limit h − = 0.That is, the top mode ε 1 is sharp at all fields in Fig. 2 in the data (first column) and in the calculation (third column), even though there are regions where it overlaps with states originating from the top continuum, as illustrated in Fig. 2, last column.
This right-most column of Fig. 2 shows the INS intensity as calculated in this section, using the same parameters as for the middle two columns but with no longitudinal mean field (H 3 = 0).Comparing Fig. 2T and P with Q and M respectively, it is seen that the above description is in remarkable qualitative agreement with the experimental results in large fields 2.5 T and 1.5 T. The agreement is especially good at 2.5 T (compare Figs. 2S  and T), which is close to the localized limit.This agreement is strongly supportive of the model presented in this section, especially given that the model is relatively simple and completely analytically tractable.At lower fields, the agreement is expected to be less good as the solitons are now more delocalized and so a perturbative treatment around the localized limit is expected to be less quantitatively accurate.Nevertheless, several key trends are still reproduced: in particular, the top mode ε 1 is captured at all fields, the middle mode ε 2 is captured down to 1 T, and the lowest mode ε − is captured down to 1.5 T.
The calculations in the right-most column of Fig. 2 include the effect of band inversion on the top two bound states, as well as the terms in H 2 , as explained in Appendix D. We note that the full Hamiltonian also contains a z magnetic field, H 3 , changing the nature of the continua.This term introduces a linear confinement, splitting the continua into confinement bound states, which cannot be captured in this calculation.Despite these important effects, the tightly bound ε bound states remain relatively unaffected by the confinement, and the qualitative predictions for the INS spectra presented in this section still hold, see Appendix C. We note that, at first order, a finite longitudinal mean field (H 3 ) would be expected to increase the energies of all the ε-bound states, which would bring the calculated dispersions in the right-most column of Fig. 2 into closer agreement with the experimental data (left-most column).

VI. CONCLUSIONS
We investigated the spectrum of the Ising chain material CoNb 2 O 6 as a function of low to intermediate transverse field in the ordered phase using inelastic neutron scattering experiments.We compared the measured spectrum to predictions based on a recently refined Hamiltonian containing all relevant sub-leading terms beyond the dominant Ising exchange and found strong quantitative agreement.We then sought a physical picture of the excitations.We found that by restricting the Hilbert space to the two-soliton subspace at first order in perturbation theory, very good agreement between the calculation and experiment was still achieved.The resulting spectrum in general has three continua and three bound states, of which only the bound states contribute significant weight to the inelastic neutron scattering intensity.In order to understand the character of the bound states, we considered the localized limit, in which the soliton hopping term on alternate bonds is zero.This occurs when the applied field matches the strength of the off-diagonal exchange.We found that the bound states in this limit are of two solitons in adjacent unit cells, stabilized by hardcore repulsion leading to a change in delocalization energy.The bound states survive well away from the localized limit, suggesting that this picture has a broader domain of validity than might initially be expected.Using this physical picture, we have been able to gain both qualitative and quantitative understanding of the low energy spectrum of CoNb  site.The diagonal matrix elements are and are consistent with the expressions given in Sec.V. Within the degenerate subspace, off-diagonal matrix elements are and Hermitian conjugate.Eigenvalues and eigenvectors are then obtained by direct diagonalization, with the dynamical correlations S xx obtained from the eigenvectors as described in Sec.V C. We term the resulting modes ε 1 and ε 2 .
We also derive the matrix elements of the next nearest neighbour terms in H 2 .Consistent with the considerations above, we only keep the diagonal matrix elements between bound states of the same type, as well as the off-diagonal matrix elements between the middle and top modes ε 0 and ε + , which show a strong mixing at the experimental parameters.The effect of the perturbation V AF , (B3), is to lower the energy of all states relative to the fully aligned (ground) state.However, as discussed in Appendix B, this energy shift is different for states with a single spin flip compared to states with at least two spin flips, corresponding to two and four satisfied antiferromagnetic bonds, respectively.This leads to matrix elements This perturbation also shifts the energies of all the continua by −2Jλ AF .The effect of the perturbation V xy AF , (B2), is to hop single spin flips by two sites.This leads to matrix elements The effect of this term on the continua is to confine the soliton pairs into a series of bound states; this effect cannot be captured within the current picture.
Finally we note that the first term in H 2 , V A , (B1), vanishes when projected to the subspace of bound states, since this term causes solitons to hop by two sites at a time.To understand the effect of this term on the continua, we note that V A , a hopping term to a neighboring unit cell, shifts the single soliton dispersion relations as ω ± −→ ω ± + Jλ A cos(kc).
In contrast to the effect of the hopping h − , V A induces the same shift in the energies of bonding / antibonding orbitals.As a result, this perturbation leads to the same energy change for the three continua, for σ, σ ′ = ±, with total and relative momenta given by k = k 1 + k 2 and q = k 1 − k 2 , respectively.Alternatively, these spectra can be obtained by applying first order perturbation theory around the localized limit, similarly to the analysis of the hopping h − presented in the main text.To this end, we first evaluate the effect of V A on the states |n, k⟩ constructed in Sec.V B. We find corresponding to an effective nearest neighbor hopping Hamiltonian for the relative coordinate n, with hopping amplitude Jλ A cos(kc/2), subject to the hard core constraint n > 0. For the top and bottom continua, the hopping amplitude due to V h− is also real, so the |q, k⟩ ± states constructed in Sec.V B are also eigenstates of V A and the change in the energies is ± ⟨q, k|V A |q, k⟩ ± ± ⟨q, k|q, k⟩ ± = 2Jλ A cos kc 2 cos qc 2 .
For the middle continuum, we consider the effect of V h− + V A on the plane wave n e iqcn/2 (|n, k⟩ 0 + + |n, k⟩ 0 − ).
We find that the perturbation corresponds to an effective nearest neighbour hopping Hamiltonian for the relative coordinate n, with complex hopping amplitude t = Jλ A cos(kc/2) + ih − sin(kc/2) = t ′ + it ′′ = |t|e iφ .This yields the dispersion The eigenstates satisfying the hard core repulsion boundary condition at n = 0 can be obtained by noting that the plane wave with relative momentum q is degenerate with the plane wave with relative momentum −q −4φ/c.Mixing these plane waves leads to the unnormalized eigenstates satisfying the hard core constraint, cell separation, n>0 e iqcn/2 − e −i(qc/2+2φ)n |n, k⟩ 0 + + |n, k⟩ 0 − .
For the plane wave n e iqcn/2 (|n, k⟩ 0 + − |n, k⟩ 0 − ), the effect of V h− + V A corresponds to an effective nearest neighbour hopping Hamiltonian for the relative coordinate n, with complex hopping amplitude t * , such that the argument above applies with φ → −φ.Thus, the effect of the perturbation V A is to add 2Jλ A cos(kc/2) cos(qc/2) to the energies of all continua, as anticipated based on the single soliton dispersion relation.We also note that the argument above yields the following two degenerate eigenstates in the presence of V h− but without V A , i.e. for φ = π/2, n>0 e icq±n/2 − e −i(cq±/2+π)n |n, k⟩ 0 + ± |n, k⟩ 0 − , with q + − q − = π.The eigenstates constructed in the main text are the symmetric / antisymmetric combinations of these eigenstates.For λ A < 0 and h − > 0 such as is found experimentally, the perturbation V A leads the top continuum to narrow and the bottom continuum to broaden, and the middle continuum to broaden around what would otherwise be the nodes.The plots in the right-most column of Fig. 2 include the effects of all terms in H 1 and H 2 , but not H 3 since it is not possible to include the effects of this last term on the continua in this framework.

J
with p L/R labeling the two site unit cells, σ L/R = 0, 1 distinguishing the even and odd sublattice and c(p L + p R )/2 being the position of the center of mass of the soliton pair.For a fixed center of mass momentum k, we obtain coupled equations for Φ (k) σ L σ R (n), defined on the half line n ≥ 0 with boundary conditions Φ = 0. We present more details on the numerical solution of these equations in Appendix A.

FIG. 5 .
FIG.5.Eigenstates of the two-soliton Hamiltonian as a function of wavevector and energy at different fields, increasing from top to bottom.In each panel, states with the character of the ε bound states have been highlighted in black.These states have been identified by being well separated from other states in energy, or, where they overlap with other states, by having high INS intensity and discarding regions with strong hybridization.Left column: Solution of H1.Three continua and three bound states can be seen at every non-zero field.ε1, ε2 and ε− are defined in Appendix D. Right column: solution of H = H1 + H2 + H3.This comparison illustrates the key effect of H3, i.e., the presence of a longitudinal mean field: all continua are split into confinement bound states, but the ε bound states are left essentially intact.This column is to be compared with the third column of Fig.2, which shows the intensities under the same conditions.
[20]6 in the low transverse field ordered phase.W. acknowledges support from a doctoral studentship funded by Lincoln College and the University of Oxford.I.L. acknowledges support from the Gordon and Betty Moore Foundation through Grant GBMF8690 to UCSB and from the National Science Foundation under Grant No. NSF PHY-1748958.D.P. acknowledges support from the Engineering and Physical Sciences Research Council grant number GR/M47249/01.L.B. was supported by the NSF CMMT program under Grant No. DMR-2116515, and by the Simons Collaboration on Ultra-Quantum Matter, which is a grant from the Simons Foundation (651440).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).The neutron scattering measurements at the ISIS Facility were supported by a beamtime allocation from the Science and Technology Facilities Council.Access to the data will be made available from Ref.[20]. ACKNOWLEDGMENTSL.