# Magnon interactions in the quantum paramagnetic phase of CoNb2O6

C. A. Gallegos  Department of Physics, University of California, Irvine,
California 92697, USA  A. L. Chernyshev  Department of Physics, University of
California, Irvine, California 92697, USA

###### Abstract

In this work, we study effects of magnon interactions in the excitation
spectrum of CoNb2O6 in the quantum paramagnetic phase in transverse field,
where the $1/S$ spin-wave theory exhibits unphysical divergences at the
critical field. We propose a self-consistent Hartree-Fock approach that
eliminates such unphysical singularities while preserving the integrity of the
singular threshold phenomena of magnon decay and spectrum renormalization that
are present in both theory and experiment. With the microscopic parameters
adopted from previous studies, this method yields a close quantitative
agreement with the available experimental data for CoNb2O6 in the relevant
regime. Insights into the general structure of the spin-anisotropic model of
CoNb2O6 and related zigzag chain materials are also provided and a discussion
of the effects of additional longitudinal field on the spectrum is given.

##  I Introduction

Quantum magnets continue to generate an enormous interest as a platform for
realizing unconventional ordered [1, 2, 3, 4, 5, 6, 7] and exotic quantum
disordered spin-liquid phases [8, 9, 10, 11], which occur due to competing
interactions between their low-energy spin degrees of freedom. The celebrated
anisotropic-exchange Kitaev model, exhibiting a spin-liquid ground state and
fractionalized excitations [12, 13, 14], has been particularly inspirational.
However, the description of real materials consistently requires exchanges
beyond the much-desired Kitaev one [15, 16, 17, 18], resulting in the models
with several bond-dependent terms, which typically favor magnetically ordered,
if unconventional, states [19, 20, 21, 22, 23, 24, 25, 26, 27, 28].

In the pursuit of the unusual physical outcomes of the bond-dependent
exchanges, recent studies have been centered on the materials with strong
spin-orbit coupling [29, 30, 31] and, specifically, on the transition-metal
compounds with Co2+ ions in an edge-sharing octahedral environment [32, 33,
34, 35, 36, 37, 38]. Cobalt niobate, CoNb2O6, is such an anisotropic-exchange
magnet, with spins forming quasi-one-dimensional ferromagnetic zigzag chains.
This material is one of the closest realizations of the Ising model, which
exhibits a paradigmatic quantum phase transition in transverse field [39]. The
field-induced transition is from the ordered phase with the domain-wall-like
excitations to the fluctuating paramagnetic phase, in which excitations are
magnon-like spin flips; see Fig. 1.

CoNb2O6 has generated further excitement by providing experimental evidence of
the bound states in its ordered phase, of the emergent E8 symmetry near the
critical field, and of the spectacular realization of the magnon decay effect
in its paramagnetic phase [40, 41, 42]. More recently, all of these phenomena
have received a consistent explanation within the microscopic spin model,
which included important bond-dependent off-diagonal exchange interactions
allowed by the crystal symmetry [43, 44, 45].

Specifically, in the paramagnetic phase of CoNb2O6, these off-diagonal
exchanges naturally yield the so-called cubic anharmonic term that couples
single-magnon excitations to the two-magnon continuum, leading to magnon
decays, the scenario confirmed by the time-dependent density-matrix
renormalization group (tDMRG) calculations [43]. The magnon decay effect is
well-documented in the isotropic and diagonal-exchange models, in which the
noncollinear states are required for the anharmonic term to occur [46, 47].
Conversely, in the presence of the off-diagonal exchanges, the anharmonic term
should be generally unavoidable even in the collinear states [48, 49, 27, 50],
which is the case of the fluctuating, nominally field polarized paramagnetic
phase of CoNb2O6 for $H\\!>\\!H_{c}$.

Figure 1: The schematic $H$–$T$ phase diagram of the 1D Ising chain in a
transverse magnetic field; $H_{c}$ is the critical field. Dashed lines
indicate the spin excitation gap.

Therefore, it is expected that the analytical insights into magnon
interactions and decays can shed further light on the important aspects of the
excitation spectrum in the paramagnetic phase of CoNb2O6 and other
anisotropic-exchange magnets. However, quantum fluctuations shift $H_{c}$ from
its classical value, leaving a wide field range inaccessible to the standard
spin-wave theory (SWT). Moreover, the $1/S$-expansion in anisotropic-exchange
models, needed to account for magnon decays, is contaminated by the unphysical
divergences at the critical field. Different methods have been proposed to
overcome similar problems in other systems, with auxiliary fields allowing for
the shifts of the phase boundaries [51, 52] and self-consistent methods
regularizing unphysical divergences [53, 50, 54].

In this work, we propose a method that allows us to explore the paramagnetic
phase of CoNb2O6 in the field range inaccessible to the standard SWT. It
naturally regularizes the unphysical $1/S$-divergences, while preserving the
integrity of the physical threshold singularities and affiliated decay
phenomena. This method combines a self-consistent Hartree-Fock approach [55,
56] with the perturbative treatment of the cubic anharmonicities. Our results
for the dynamical structure factor using microscopic parameters suggested in
Ref. [44] agree closely with the inelastic neutron scattering data [42, 43]
that were previously reproduced by the tDMRG [43]. Furthermore, we investigate
the effects of the additional longitudinal fields in the spectrum of CoNb2O6
in the paramagnetic phase and make predictions of magnon band gaps and
associated Van Hove singularities. We also provide insights into the general
structure of the spin-anisotropic model for this and related materials.

The paper is organized as follows. In Section II, we introduce the spin
Hamiltonian for CoNb2O6, relate it to the other anisotropic models, and
discuss phenomenological constraints. In Section III, we present the standard
$1/S$ spin-wave expansion, demonstrate the unphysical divergences in it, and
describe the self-consistent method that regularizes them. In Section IV, we
compare our results with the experimental data. In Section V, we present our
predictions of the effects of additional longitudinal fields in the spectrum
of CoNb2O6. We conclude by summarizing our results in Section VI. Appendixes
provide technical details.

##  II Model

In this Section, we introduce the anisotropic-exchange model that should
describe magnetic properties of CoNb2O6 and related quasi-1D materials.
Following the prior analysis [44, 43], we use the space group of the crystal
structure, provide a connection of this model to the broader class of
anisotropic-exchange models, and discuss phenomenological constraints on the
spin Hamiltonian given in Refs. [44, 43].

###  II.1 Crystal structure

The crystal structure of CoNb2O6 is orthorhombic, space group $Pbcn$. The
combination of the crystal-field and spin-orbit coupling splits the
$j\\!=\\!3/2$ multiplet of the Co2+, leading to an effective spin-$1/2$ ground
state on each magnetic site [57]. The magnetic Co2+ ions are arranged in 1D
zigzag chains oriented along the crystallographic $c$ axis with the staggered
displacement along the $b$ axis; see Fig. 2(a). In the basal $ab$ plane, Co2+
spins form a weakly-coupled deformed triangular lattice; see Ref. [41] for
details. The schematic representation of the isolated spin-$1/2$ zigzag chain
is shown in Fig. 2(a) together with the crystallographic $\\{a,b,c\\}$ axes.

At low temperatures, Co2+ moments in each chain order ferromagnetically,
pointing along the Ising easy-axis, which lies in the $ac$ plane at an angle
$\gamma\\!\approx\\!30^{\circ}$ to the $c$ axis [57]. Therefore, it is natural
to introduce another reference frame, referred to as the laboratory frame
$\\{x_{0},y_{0},z_{0}\\}$, obtained by a rotation of the crystallographic
frame about the $b$ axis by the angle $\gamma$, so that $z_{0}$ is aligned
with the Ising direction; see Fig. 2(b).

Figure 2: (a) The schematic representation of a segment of the zigzag spin
chain, with the crystallographic $\\{a,b,c\\}$ axes, bond-dependent angles
$\varphi_{\alpha}$, nearest-neighbor vectors $\bm{\delta}_{\alpha}$, lattice
constant ${\sf c}$, width of the chain ${\sf b}$, glide $ac$ plane, and $bc$
plane of the zigzag structure. Dashed lines in $b$ direction are the imaginary
missing bonds of the hypothetical honeycomb lattice, see text. (b) The
laboratory reference frame $\\{x_{0},y_{0},z_{0}\\}$ with $z_{0}$ axis along
the zero-field Ising direction of spins and angle $\gamma$ in the $ac$ plane.
(c) The “local” reference frame $\\{x,y,z\\}$ in the paramagnetic phase with
$z$ axis along the transverse-field-induced spin orientation for
$H\\!>\\!H_{c}$.

###  II.2 Symmetries and the nearest-neighbor model

Here, we consider the 1D zigzag spin-chain model with interactions only
between the nearest-neighbor sites. Given the translational invariance, the
most general nearest-neighbor spin Hamiltonian of such a chain is

$\displaystyle\mathcal{\hat{H}}_{1}$ $\displaystyle\\!=\\!\sum_{\langle
ij\rangle}\textbf{S}^{\text{T}}_{i}\hat{\bm{J}}_{\alpha}\textbf{S}_{j},$ (1)

where $\textbf{S}_{i}^{\text{T}}\\!=\\!(S_{i}^{x},S_{i}^{y},S_{i}^{z})$,
$\langle ij\rangle$ denotes summation over the nearest-neighbor bonds,
$\alpha\\!=\\!1,2$ numerates the two distinct bonds of the zigzag structure
with the nearest-neighbor vectors $\bm{\delta}_{1(2)}$ depicted in Fig. 2(a),
and $\hat{\bm{J}}_{1(2)}$ being their respective $3\\!\times\\!3$ exchange
matrices. At this stage, the two exchange matrices in the model (1) have
eighteen independent parameters in total.

The number of independent parameters in the model (1) is reduced by the space
group symmetry of the lattice. The effect of these symmetries on the form of
the exchange matrices $\hat{\bm{J}}_{\alpha}$ have been thoroughly discussed
in Ref. [43]. Here, we provide a complementary derivation.

CoNb2O6 has two space-group symmetries, the bond-center inversion of the
nearest-neighbor bond and the glide symmetry. The inversion with respect to
the bond center transposes individual exchange matrices
$\hat{\bm{J}}_{\alpha}$, but must leave them invariant, permitting only
symmetric off-diagonal terms. This reduces the number of independent
parameters in the model (1) to six per bond.

The glide symmetry consists of the spatial reflection in the $ac$ plane,
followed by a translation by half of the unit cell ${\sf c}_{0}\\!=\\!{\sf
c}/2$; see Fig. 2(a). The spatial reflection flips the sign of the spin
components that are parallel to the $ac$ plane, leaving the $b$-component
intact. The half-translation completes the space-group operation, but swaps
$\hat{\bm{J}}_{1}$ and $\hat{\bm{J}}_{2}$, yielding the exchange matrices for
the two bonds in the crystallographic $\\{a,b,c\\}$ reference frame given by

$\hat{\bm{J}}_{\alpha}\\!=\\!\begin{pmatrix}{J}_{aa}&(-1)^{\alpha}{J}_{ab}&{J}_{ac}\\\
(-1)^{\alpha}{J}_{ab}&{J}_{bb}&(-1)^{\alpha}{J}_{bc}\\\
{J}_{ac}&(-1)^{\alpha}{J}_{bc}&{J}_{cc}\\\ \end{pmatrix}.$ (2)

Thus, the nearest-neighbor model (1) has only six independent parameters,
$\\{J_{aa},J_{bb},J_{cc},J_{ab},J_{ac},J_{bc}\\}$.

An important aspect of the exchange matrices in (2) is the presence of the two
off-diagonal staggered terms that alternate between the zigzag bonds. Such a
bond-dependence suggests a broad relation of the model for CoNb2O6 with the
other well-known forms of the anisotropic-exchange models, discussed next.

###  II.3 Alternative parameterizations

Given the bond-dependence of the exchange matrices in (2), it is tempting to
establish a connection between the zigzag chain model and the much-studied
bond-dependent models on the honeycomb and other lattices. To make this
parallel more explicit geometrically, one can perceive the 1D zigzag structure
as an element of a hypothetical honeycomb lattice, which is missing bonds in
one direction [58, 59, 60]. For CoNb2O6, one can introduce the imaginary
missing bonds in the $b$ direction; see Fig. 2(a). As is noted below, the
mutual 2D arrangement of the chains in the $bc$ plane of CoNb2O6 does not
correspond to the honeycomb lattice, but an important symmetry that is needed
to make such a construct possible is present. However, there are a few nuances
in the discussed connection that are worth highlighting.

First, the angles of the nearest-neighbor vectors $\bm{\delta}_{\alpha}$ with
the imaginary missing bonds shown in Fig. 2(a), $\varphi_{1,2}\\!=\\!\mp
127{\degree}$, are close but not equal to those in an ideal honeycomb lattice.
Second, unlike in the honeycomb lattice, the physical bonds of the zigzag
chain are not the $C_{2}$-symmetry axes, or, alternatively, the zigzag chain
has only one of the three glide planes of the honeycomb lattice. While the
imaginary bonds are the true $C_{2}$-symmetry axes, the $\pi$-rotation in them
is equivalent to a combination of the glide and bond-center inversion
symmetries discussed above, providing no further restrictions on the
parameters of the model in Eq. (2).

Curiously, the true 2D arrangement of the chains in the $bc$ plane of CoNb2O6
is that of a distorted centered rectangular lattice, see Fig. 3, which has the
$C_{2}$-symmetry axis for the imaginary bonds.

Figure 3: The crystal structure of CoNb2O6 as seen in the $bc$ plane, (a) with
and (b) without the oxygen octahedral environment. The Nb ions are not shown.

With these insights, the model (2) can be straightforwardly cast into the
“ice-like” form [3]. Within this parametrization, one remains in the same
crystallographic reference frame $\\{a,b,c\\}$, but the diagonal elements in
(2) are rewritten as

$J_{bb}+J_{cc}\\!=\\!2J,\
J_{bb}-J_{cc}\\!=\\!4J_{\pm\pm}\cos{\varphi_{\alpha}},\ J_{aa}\\!=\\!\Delta
J,$ (3)

where $\Delta$ is the $XXZ$ anisotropy parameter, with $a$ being anisotropy
axis, and $\varphi_{\alpha}$ are the bond angles in Fig. 2(a). The off-
diagonal terms can be rewritten as
$J_{ac}\\!=\\!J_{z\pm}\cos{\varphi_{\alpha}}$ and

$\displaystyle(-1)^{\alpha}\\!J_{ab}\\!=\\!-\lambda_{z}J_{z\pm}\\!\sin{\varphi_{\alpha}},\,(-1)^{\alpha}\\!J_{bc}\\!=\\!2\lambda_{\pm}J_{\pm\pm}\\!\sin{\varphi_{\alpha}},$
(4)

thus encoding the staggered nature of the bond-dependent terms in that of the
bond angles, $\varphi_{1}\\!=\\!-\varphi_{2}$; see Fig. 2(a). This converts
the bond-dependent exchange matrix in (2) to

$\hat{\bm{J}}_{\alpha}\\!=\\!\begin{pmatrix}\Delta
J&-\lambda_{z}J_{z\pm}s_{\alpha}&J_{z\pm}c_{\alpha}\\\
-\lambda_{z}J_{z\pm}s_{\alpha}&J+2J_{\pm\pm}c_{\alpha}&2\lambda_{\pm}J_{\pm\pm}s_{\alpha}\\\
J_{z\pm}c_{\alpha}&2\lambda_{\pm}J_{\pm\pm}s_{\alpha}&J-2J_{\pm\pm}c_{\alpha}\\\
\end{pmatrix},$ (5)

where the notations $c_{\alpha}\\!\equiv\\!\cos{\varphi}_{\alpha}$ and
$s_{\alpha}\\!\equiv\\!\sin{\varphi}_{\alpha}$ are used for brevity. The form
(5) provides an alternative parameterization to the exchange matrix,
translating the set of six independent variables in (2) to
$\\{J,\Delta,J_{\pm\pm},J_{z\pm},\lambda_{\pm},\lambda_{z}\\}$, see Appendix
A.1.

This form in (5) is similar to the anisotropic-exchange matrices for the
triangular and honeycomb lattices within the same “ice-like” parametrization,
up to a cyclic permutation of the axes [26, 50], but it is enriched by the two
additional independent terms, which are introduced as the multiplicative
factors, $\lambda_{\pm}$ and $\lambda_{z}$. The presence of these extra terms
is due to the lower symmetry of the zigzag chains discussed above. An obvious
utility of the form (5) is that one can straightforwardly characterize the
deviation of the zigzag model from the more symmetric honeycomb-lattice case
by using the actual parameters proposed for CoNb2O6 and examining the
differences of $\lambda_{\pm}$ and $\lambda_{z}$ from unity.

We also note that recently, an attempt to introduce the bond-dependent
exchanges of the Kitaev honeycomb-lattice model to describe CoNb2O6 was made
in the form of the twisted Kitaev-chain Hamiltonian; see Ref. [61]. This
Hamiltonian corresponds to an interpolation between the Ising chain and the 1D
Kitaev model [62]. However, a limited number of independent exchange terms in
this model restricts its ability to describe quantitatively generic
anisotropic-exchange zigzag chain materials and, specifically, CoNb2O6 [44].
Naturally, a more complete description should be achievable within this
Kitaev-like parametrization, but it would require an extended version of the
model with all symmetry-allowed terms present in the exchange matrices,
corresponding to a generalized Kitaev-Heisenberg chain [58, 59, 60].

Next, we discuss phenomenological constraints on the model (2) for CoNb2O6.

###  II.4 Phenomenological constraint

In principle, having taken advantage of all of the lattice symmetries, the
full six-parameter nearest-neighbor model (1) with the exchange matrices from
(2) and minimal additional further neighbor and three-dimensional terms should
be used to provide the best fit of the experimental data in order to determine
the actual values of these parameters for a specific material.

However, in the most comprehensive studies of CoNb2O6 in Refs. [43, 44], the
number of independent terms in the nearest-neighbor exchange matrices (2) has
been reduced to four by utilizing a phenomenological constraint before the
parameter fitting procedure.

####  II.4.1 Ising axis direction

Since the zero-field magnetic order in CoNb2O6 has a preferred direction, it
is natural to rotate the crystallographic reference frame to align one of the
axes ($z_{0}$) with the observed Ising axis of the spins. It is done by a
rotation of the $\\{a,b,c\\}$ axes by $\gamma$ about the $b$ axis; see Fig.
2(b). Then, the exchange matrices (2) are transformed to this laboratory
reference frame $\\{x_{0},y_{0},z_{0}\\}$ via
$\hat{{\textbf{J}}}_{\alpha}\\!=\\!\hat{\textbf{R}}_{\gamma}\hat{{\bm{J}}}_{\alpha}\hat{\textbf{R}}_{\gamma}^{\text{T}}$,
where $\hat{\textbf{R}}_{\gamma}$ is the rotation matrix and

$\hat{\textbf{J}}_{\alpha}\\!=\\!\begin{pmatrix}J_{x_{0}x_{0}}&(-1)^{\alpha}J_{x_{0}y_{0}}&J_{x_{0}z_{0}}\\\
(-1)^{\alpha}J_{x_{0}y_{0}}&J_{y_{0}y_{0}}&(-1)^{\alpha}J_{y_{0}z_{0}}\\\
J_{x_{0}z_{0}}&(-1)^{\alpha}J_{y_{0}z_{0}}&{J}_{z_{0}z_{0}}\\\ \end{pmatrix},$
(6)

with the relations of the $\\{x_{0},y_{0},z_{0}\\}$ exchanges to the ones in
the $\\{a,b,c\\}$ frame given in Appendix A.1. One can notice that matrices in
(6) retain the structure of (2).

Since the Ising direction is minimizing the energy of the zero-field spin
configuration, it provides an implicit phenomenological constraint on the
matrix elements of $\hat{\textbf{J}}_{\alpha}$ in (6), such that the spins in
the ground state of the model should stay aligned along the $z_{0}$ axis.

The essence of the approach proposed in Refs. [43, 44] is to impose such a
constraint explicitly by eliminating all individual terms in (6) that generate
an unphysical tilt of spins away from the $z_{0}$ axis. One of such offending
terms is $J_{x_{0}z_{0}}$. Since it creates an $x_{0}$-tilt of spins in the
$x_{0}z_{0}$ plane already in the classical limit of the model, it is rendered
zero in this approach. Curiously, the $J_{y_{0}z_{0}}$-term does not provide a
$y_{0}$-tilt because of its staggered nature stemming from the glide symmetry
of the lattice.

Less obviously, in the quantum case, the $x_{0}$-tilt is also generated by the
combination of the two staggered terms, $J_{x_{0}y_{0}}$ and $J_{y_{0}z_{0}}$,
as we demonstrate below. The $J_{y_{0}z_{0}}$-term was found crucial for the
CoNb2O6 phenomenology as the key microscopic source of the domain-wall
dispersion observed in the ordered phase [43, 44]. Then, it follows that the
only way to eliminate the unphysical $x_{0}$-tilt completely is to vanish
$J_{x_{0}y_{0}}$-term, yielding the four-parameter exchange matrix advocated
in Refs. [43, 44],

$\hat{\textbf{J}}_{\alpha}\\!=\\!\begin{pmatrix}J_{x_{0}x_{0}}&0&0\\\
0&J_{y_{0}y_{0}}&(-1)^{\alpha}J_{y_{0}z_{0}}\\\
0&(-1)^{\alpha}J_{y_{0}z_{0}}&{J}_{z_{0}z_{0}}\\\ \end{pmatrix}.$ (7)

While we will adopt this form of the exchange matrix in the main part of the
present work below, the following note is in order. Although the approach of
Refs. [43, 44] is simple, seemingly unambiguous, and potentially generic, it
is not without a caveat.

One may suspect that such an approach is overconstraining, because a single
phenomenological constraint is used to eliminate two symmetry-allowed terms
from the exchange matrix. Instead, both offending terms, $J_{x_{0}z_{0}}$ and
$J_{x_{0}y_{0}}$, may be allowed to be non-zero, but exactly compensating each
others’ spin tilting and leaving the physical Ising direction intact. Of
course, the technical implementation of such an indirect constraint as a part
of the parameter-fitting procedure is more challenging, so the approach of
making $J_{x_{0}z_{0}}\\!=\\!J_{x_{0}y_{0}}\\!=\\!0$ can be taken as a mild
assumption in the search of a minimal model.

To vindicate the assumption of Refs. [43, 44] in the case of CoNb2O6 further,
we note that the compensating tilts from $J_{x_{0}z_{0}}$ and
$J_{x_{0}y_{0}}(J_{y_{0}z_{0}})$ terms appear in different orders of the
quasiclassical theory. $J_{x_{0}z_{0}}$ creates a tilt already in the
classical limit of the model, while the tilt due to $J_{x_{0}y_{0}}$ is a
purely quantum effect. Given this hierarchy and using the fact that the off-
diagonal exchanges in CoNb2O6 are secondary to the main Ising term, below we
provide a perturbative consideration of the effects of the “residual”
$J_{x_{0}z_{0}}$ and $J_{x_{0}y_{0}}$ terms.

####  II.4.2 Perturbative consideration

Coming back momentarily to the exchange matrix in the crystallographic
reference frame (2), a straightforward minimization of the classical energy of
the model yields the tilt angle of spins away from the $c$ axis as

$\tan
2\widetilde{\gamma}=\frac{2J_{ac}}{J_{cc}-J_{aa}}=-\frac{2J_{z\pm}\cos\varphi_{\alpha}}{(\Delta-1)J+2J_{\pm\pm}\cos\varphi_{\alpha}},$
(8)

given here for both parametrizations of the exchange matrix in (2) and (5).
The latter illustrates one of the broader perspectives provided by the form in
(5) as the $J_{z\pm}$ term is known to produce such an out-of-the-plane tilt
in the previously discussed models [63, 26, 50].

Then, within the classical approximation, one would equate the tilt angle
$\widetilde{\gamma}$ to its experimentally observed value $\gamma$, thus using
the preferred direction of the magnetic order in CoNb2O6 shown in Fig. 2(b) as
a phenomenological constraint that provides a relation between exchanges given
by Eq. (8). As a result, the number of independent terms in the nearest-
neighbor exchange matrix would be reduced to five, setting
$J_{x_{0}z_{0}}\\!\equiv\\!0$ in (6).

However, quantum fluctuations can renormalize the tilt angle of the ordered
magnetic moment, producing deviations from the classical result (8). In other
words, if one would calculate the angle between the Ising $z_{0}$ and $c$ axes
in the quantum case with $\hat{{\textbf{J}}}_{\alpha}$ in (6) and
$J_{x_{0}z_{0}}\\!=\\!0$, it would generally deviate from $\gamma$.

This quantum renormalization can be accessed perturbatively by considering
virtual spin-flip processes that are generated by the staggered terms
$J_{x_{0}y_{0}}$ and $J_{y_{0}z_{0}}$. Using the real-space perturbation
theory [64, 65, 66] for the $S\\!=\\!1/2$ model in (6) with only the main
Ising and staggered terms, we derive the tilt angle in the second order of the
theory as given by

$\delta\gamma\approx-\frac{J_{x_{0}y_{0}}J_{y_{0}z_{0}}}{J_{z_{0}z_{0}}^{2}},$
(9)

which is supported by our DMRG calculations [67], with the details for both
deferred to Appendix A.2.

We can further assert that the higher-order corrections to the tilt angle also
require both staggered terms, because they have to cancel their symmetry-
related staggered form as in (9). Moreover, the higher-order corrections to
(9) need to carry odd powers of each of the staggered terms because they
generate different number of spin flips, as can also be verified numerically;
see Appendix A.2.

Thus, from the quasiclassical perspective, the choice of
$J_{x_{0}y_{0}}\\!=\\!0$ made in Refs. [43, 44] automatically renders all
quantum corrections to the classical Ising axis angle $\widetilde{\gamma}$
equal to zero, leaving it equal to the experimental value $\gamma$ by
construction. For the parametrizations of the exchange matrices in (2) and
(5), the choice $J_{x_{0}y_{0}}\\!=\\!0$ provides another relation between the
components of the exchange matrices that reads

$\tan\gamma=\frac{J_{ab}}{J_{bc}}=-\frac{\lambda_{z}J_{z\pm}}{2\lambda_{\pm}J_{\pm\pm}}.$
(10)

Altogether, for the nearest-neighbor model written in the laboratory reference
frame $\\{x_{0},y_{0},z_{0}\\}$, this results in the four-parameter exchange
matrix given in (7).

Lastly, one can use the perturbative consideration for the tilt angle in (9)
together with the numerically precise DMRG calculations in order to quantify
the potential values of the “residual,” mutually compensating $J_{x_{0}z_{0}}$
and $J_{x_{0}y_{0}}$ terms in the quantum $S\\!=\\!1/2$ model of CoNb2O6, if
these terms are allowed to deviate from zero. A straightforward derivation
gives the relation between such terms that would leave the Ising axis
direction intact,

$J_{x_{0}z_{0}}=J_{x_{0}y_{0}}\,\left(\frac{J_{y_{0}z_{0}}}{J_{z_{0}z_{0}}}\right),$
(11)

which explicates the different order of their corresponding effects on the
spin orientation in the quasiclassical expansion. We note that this result is
obtained for a simplified model as is Eq. (9).

In Fig. 4, we show this dependence for the choice of $J_{y_{0}z_{0}}$ and
$J_{z_{0}z_{0}}$ values that correspond to CoNb2O6; see Sec. II.5 below. It is
plotted together with the DMRG results for the full model using the best-fit
parameters discussed in the next Section. According to Ref. [43], the
$J_{x_{0}y_{0}}$ term in CoNb2O6 should be small as it produces the momentum
space periodicity of the domain-wall excitations in the ordered phase that is
different from the observed one. As is shown in Fig. 4, this should render the
limit on the residual $J_{x_{0}z_{0}}$ to nearly zero, thus providing a
further partial exoneration to the approach of Refs. [43, 44].

Figure 4: $J_{x_{0}z_{0}}$ vs $J_{x_{0}y_{0}}$ for the CoNb2O6 parameters in
Sec. II.5. Line is the real-space perturbation theory (RSPT) result (11) for
the simplified model with only the $J_{z_{0}z_{0}}$, $J_{x_{0}y_{0}}$, and
$J_{y_{0}z_{0}}$ terms; symbols are DMRG results for the full model, see
Appendix A.2 for details.

###  II.5 Best-fit parameters

An extensive comparison of the experimental data with the numerical modeling
carried out in Ref. [44] has resulted in the set of the best-fit parameters
for CoNb2O6. Translating them to the notations of our Eq. (7), the nearest-
neighbor exchanges are

$\displaystyle J_{x_{0}x_{0}}\\!$ $\displaystyle=\\!-0.57(2)\ \text{meV},\
J_{y_{0}y_{0}}\\!=\\!-0.67(2)\ \text{meV},$ $\displaystyle J_{z_{0}z_{0}}\\!$
$\displaystyle=\\!-2.48(2)\ \text{meV},\ J_{y_{0}z_{0}}\\!=\\!-0.56(1)\
\text{meV}.$ (12)

We note that in Refs. [43, 44] a different parametrization has been used for
the two diagonal exchanges:
$J_{x_{0}x_{0}}\\!=\\!J_{z_{0}z_{0}}(\lambda_{s}\\!+\\!\lambda_{a})$ and
$J_{y_{0}y_{0}}\\!=\\!J_{z_{0}z_{0}}(\lambda_{s}\\!-\\!\lambda_{a})$.

For completeness, we also translate these numerical values to that of the
exchange parameters in the original crystallographic $\\{a,b,c\\}$ frame in
Eq. (2) and to the “ice-like” parametrization of Eq. (5); see Appendix A.1.
For the “ice-like” form (5), the best-fit parameters
$\\{J,\Delta,J_{\pm\pm},J_{z\pm},\lambda_{\pm},\lambda_{z}\\}$ are given by

$\displaystyle J\\!$ $\displaystyle=\\!-1.34(1)\ \text{meV},\
\Delta\\!=\\!0.78(1),\ J_{\pm\pm}\\!=\\!-0.55(1)\ \text{meV},$ $\displaystyle
J_{z\pm}\\!$ $\displaystyle=\\!1.37(2)\ \text{meV},\
\lambda_{z}\\!=\\!0.26(1),\ \lambda_{\pm}\\!=\\!0.55(1),$ (13)

where the last two parameters, $\lambda_{z}$ and $\lambda_{\pm}$, quantify the
substantial degree to which the zigzag chain differs from the hypothetical
honeycomb lattice, where $\lambda_{z}\\!=\\!\lambda_{\pm}\\!=\\!1$. Clearly,
the diagonal $XXZ$ exchanges are ferromagnetic, and the values of $J_{\pm\pm}$
and $J_{z\pm}$ underscore the pronounced anisotropic nature of CoNb2O6.
Interestingly, the value of $\Delta\\!<\\!1$ implies that in this
parametrization, CoNb2O6 can be regarded as an easy-plane anisotropic-exchange
magnet, thus establishing a connection to the other members of the cobaltate
family [68, 36, 37]. Moreover, the dominant $J_{z\pm}$ and smaller
$J_{\pm\pm}$ anisotropies are reminiscent of the other transition-metal
anisotropic-exchange materials, such as $\alpha$-RuCl3 [27].

As is stated in Ref. [44], the nearest-neighbor model (1) with the exchange
matrix in (7), needs to be supplemented by the next-nearest-neighbor $XXZ$
term

$\hat{\mathcal{H}}_{2}\\!=\\!\sum_{i}\Big{\\{}J_{2}\left(S_{i}^{x_{0}}S_{i+2}^{x_{0}}+S_{i}^{y_{0}}S_{i+2}^{y_{0}}\right)\\!+\\!J_{2z_{0}}S_{i}^{z_{0}}S_{i+2}^{z_{0}}\Big{\\}},$
(14)

providing a consistent quantitative agreement with the experimental data for
the excitation spectrum of CoNb2O6 in different regions of its phase diagram
and for the fields applied in the transverse and longitudinal directions. The
best parameter choice for the next-nearest-neighbor part of the Hamiltonian
(14) is

$\displaystyle J_{2}\\!$ $\displaystyle=\\!0.077(3)\ \text{meV},\
J_{2z_{0}}\\!=\\!0.19(1)\ \text{meV},$ $\displaystyle g_{x_{0}}\\!$
$\displaystyle=\\!3.29(6),\ g_{y_{0}}\\!=\\!3.32(2),\
g_{z_{0}}\\!=\\!6.90(5),$ (15)

where we have also listed the principal moments of the g-tensor. Below, we
will use the best-fit parameter values in (12) and (15) to compare our results
with the neutron scattering data in the paramagnetic phase of CoNb2O6 [42,
43]. We will also follow prior works [44] in a simplifying assumption that the
g-tensor is diagonal in the laboratory frame.

We also point out that one may need to include small interchain couplings if
considering the full 3D spectrum, or use an effective longitudinal field to
account for their confining effect in the ordered phase [41, 43, 44, 45].
However, in this work, we focus on the spectrum in the field-induced
paramagnetic phase of CoNb2O6, in which spins are aligned in the transverse
($b\\!=\\!y_{0}$) direction, making the effect of the interchain couplings
negligible.

##  III Self-consistent spin-wave theory

In this Section, we briefly outline the steps of the $1/S$ spin-wave expansion
as applied to the paramagnetic phase of CoNb2O6, demonstrate the problem of
the unphysical divergences in it, and describe the self-consistent Hartree-
Fock method that regularizes them.

###  III.1 Critical field and Hamiltonian in local axes

The paramagnetic phase in CoNb2O6 is induced by a transverse magnetic field.
To model it, the six-parameter spin Hamiltonian
$\hat{\mathcal{H}}_{1}+\hat{\mathcal{H}}_{2}$ from Eqs. (1), (7), and (14) has
to be augmented by the transverse-field term

$\hat{\mathcal{H}}_{\perp}=g_{y_{0}}\mu_{B}B\sum_{i}S_{i}^{y_{0}},$ (16)

with the field $B$ along the high-symmetry $y_{0}(b)$ axis. Using the
classical energy consideration detailed in Appendix B with exchanges and
relevant g-tensor component from Eqs. (12) and (15) gives the classical
critical field for the transition to the paramagnetic phase

$H_{c}=2S(J_{y_{0}y_{0}}\\!-\\!J_{z_{0}z_{0}}\\!+\\!J_{2}\\!-\\!J_{2z_{0}}),\
B_{c}\approx 8.8(1)\ \mathrm{T},$ (17)

where the field in the energy units, $H\\!=\\!g_{y_{0}}\mu_{B}B$, is
introduced. For the sake of the future discussion, we note that this critical
field is considerably larger than the one found by the DMRG in the 1D model,
$B_{c}^{\text{1D}}\\!\approx\\!4.5\ \mathrm{T}$, see Section IV, suggesting
strong renormalization due to quantum effects.

For the spin-wave expansion in the paramagnetic phase, we perform a rotation
from the laboratory $\\{x_{0},y_{0},z_{0}\\}$ to the local reference frame
$\\{x,y,z\\}$, depicted in Fig. 2(c), aligning the local quantization axis $z$
with the direction of the field. This leads to the cyclic permutation of the
spin components

$\left(S_{i}^{x_{0}},S_{i}^{y_{0}},S_{i}^{z_{0}}\right)_{\text{lab}}=\left(S_{i}^{y},S_{i}^{z},S_{i}^{x}\right)_{\text{loc}}.$
(18)

After the transformation (18), it is convenient to divide the Hamiltonian into
two parts, referred to as the even and the odd, in order to separate even and
odd powers of the bosonic operators in the subsequent spin bosonization. Using
the Hamiltonian in Eqs. (1), (7), (14), and (16), the even term reads

$\displaystyle\hat{\mathcal{H}}_{\text{even}}$
$\displaystyle\\!=\\!\sum_{i}\Big{\\{}J_{z_{0}z_{0}}S_{i}^{x}S_{i+1}^{x}+J_{x_{0}x_{0}}S_{i}^{y}S_{i+1}^{y}+J_{y_{0}y_{0}}S_{i}^{z}S_{i+1}^{z}$
$\displaystyle+J_{2z_{0}}S_{i}^{x}S_{i+2}^{x}+J_{2}\left(S_{i}^{y}S_{i+2}^{y}+S_{i}^{z}S_{i+2}^{z}\right)\\!-\\!HS_{i}^{z}\Big{\\}},$
(19)

while the odd part is given by

$\displaystyle\hat{\mathcal{H}}_{\text{odd}}=J_{y_{0}z_{0}}\sum_{i}e^{i{\bf
Q}{\bf r}_{i}}\left(S_{i}^{x}S_{i+1}^{z}+S_{i}^{z}S_{i+1}^{x}\right).$ (20)

In the latter, the factor $e^{i{\bf Q}{\bf r}_{i}}\\!=\\!(-1)^{i}$ replicates
the staggered structure of the $J_{y_{0}z_{0}}$-term, where ${\bf
Q}\\!=\\!2\pi\hat{\textbf{{\sf c}}}/{\sf c}$ is the reciprocal lattice vector
of the zigzag chain, ${\sf c}$ is its lattice constant, and $\hat{\textbf{{\sf
c}}}$ is the unit vector along the $c$ axis in Fig. 2(a). As is discussed
below, the relevant unit cell is smaller, with the lattice constant ${\sf
c}_{0}\\!=\\!{\sf c}/2$ and the reciprocal lattice vector ${\bf G}\\!=\\!2{\bf
Q}$.

Clearly, in the absence of the odd part (20), there would be no memory of the
zigzag structure left in the spin model, as the even part (19) describes a
“simple” Ising-like chain with the transverse field term and second-neighbor
exchanges, but no bond-dependent terms. Since (19) yields the linear spin-wave
theory, one can anticipate that it will have only a single bosonic branch,
with no zone-folding in the reciprocal space from the zigzag structure of
CoNb2O6.

On the other hand, the odd part of the model (20) arises precisely from such
bond-dependent terms. However, in the paramagnetic phase, it contributes only
to the nonlinear, anharmonic coupling of the spin flips, bringing about an
important ${\bf Q}$-shift of the two-magnon continuum that couples to the
single-magnon branch. This feature of the spin model of the zigzag chains in
CoNb2O6 has been recognized and thoroughly discussed in Ref. [43] as crucial
for explaining puzzling kinematics of the observed magnon decays.

We also note that a similar structure of the theory was recently discussed in
Ref. [50] in the context of the easy-plane honeycomb-lattice model with bond-
dependent exchanges, underscoring the connection of the present consideration
to a broader class of models and materials with spin-orbit-generated
anisotropic exchanges.

###  III.2 Linear spin-wave theory

The harmonic, or linear spin-wave theory (LSWT) order of the $1/S$-expansion
about the classical ground state is obtained via the standard Holstein-
Primakoff (HP) bosonization of spin operators in the local reference frame:
$S_{i}^{z}\\!=\\!S-n_{i}$ and, to the lowest order,
$S_{i}^{+}\\!\approx\\!\sqrt{2S}a_{i}$.

In the field-polarized paramagnetic phase of CoNb2O6 considered here, it is
only the even part of the Hamiltonian (19) that contributes to LSWT. As is
discussed above, this part of the Hamiltonian is invariant to the translations
by ${\sf c}_{0}\\!=\\!{\sf c}/2$, that is, half of the primitive lattice
vector of the zigzag chain. In other words, spin states on all sites of the
chain are equivalent and only one bosonic species needs to be introduced.
Using the HP bosonization in (19) and the standard Fourier transformation

$a_{i}=\frac{1}{\sqrt{N}}\sum_{\bf k}e^{i{\bf k}{\bf r}_{i}}a_{\bf k},$ (21)

where $N$ is the number of lattice sites in the chain, we obtain the LSWT
Hamiltonian

$\mathcal{\hat{H}}^{(2)}=\sum_{\bf k}\Big{\\{}A_{\bf k}a^{\dagger}_{\bf
k}a^{\phantom{{\dagger}}}_{\bf k}-\frac{B_{\bf k}}{2}\left(a^{\dagger}_{{\bf
k}}a^{\dagger}_{-{\bf k}}+\text{H.c.}\right)\Big{\\}},$ (22)

where $A_{\bf k}$ and $B_{\bf k}$ are

$\displaystyle A_{\bf k}=H-2S(J_{y_{0}y_{0}}+J_{2})+SJ_{{\bf k}+},\ B_{\bf
k}=SJ_{{\bf k}-},$ (23)

with

$\displaystyle J_{{\bf k}\pm}=\left(J_{x_{0}x_{0}}\pm
J_{z_{0}z_{0}}\right)\gamma_{\bf k}^{(1)}+\left(J_{2}\pm
J_{2z_{0}}\right)\gamma_{\bf k}^{(2)},$ (24)

and the nearest- and next-nearest-neighbor hopping amplitudes $\gamma_{\bf
k}^{(n)}\\!=\\!\cos(nk_{c}{\sf c}_{0})$, where $n\\!=\\!1(2)$,
$k_{c}\\!=\\!{\bf k}\hat{\textbf{{\sf c}}}$ is a projection of the momentum
${\bf k}$ on the chain direction, and ${\sf c}_{0}\\!=\\!{\sf c}/2$, with
${\sf c}$ being the lattice constant of the zigzag chain, as before; see Fig.
2(a) and Appendix B for more details and Ref. [41] for similar expressions.

The LSWT Hamiltonian (22) is diagonalized by a textbook Bogolyubov
transformation, $a_{\bf k}\\!=\\!u_{\bf k}b_{\bf
k}^{\phantom{{\dagger}}}\\!+v_{\bf k}b_{-{\bf k}}^{\dagger}$, with $u_{\bf
k}^{2}+v_{\bf k}^{2}\\!=\\!A_{\bf k}/\varepsilon_{\bf k}$, $2u_{\bf k}v_{\bf
k}\\!=\\!B_{\bf k}/\varepsilon_{\bf k}$, and magnon energy

$\varepsilon_{\bf k}=\sqrt{A_{\bf k}^{2}-B_{\bf k}^{2}}.$ (25)

The excitation gap of the LSWT spectrum (25) at the $\bm{\Gamma}$ point $({\bf
k}\\!=\\!0)$ is given by

$\Delta_{0}=\sqrt{(H-H_{c})(H+2S(J_{x_{0}x_{0}}-J_{y_{0}y_{0}}))},$ (26)

vanishing at the critical field $H_{c}$ in (17), as expected. Given the low
spin-symmetry of the model, the spectrum in (25) has a relativistic form near
$\bm{\Gamma}$, with $\varepsilon_{\bf k}\\!\propto\\!|{\bf k}|$ at
$H\\!=\\!H_{c}$, the behavior that will be important for the unphysical
divergences discussed below.

Figure 5: (a) Decay, (b) source, (c) Hartree-Fock, and (d) self-consistent
Hartree-Fock self-energies.

###  III.3 Non-linear spin-wave theory and divergences

The $1/S$-expansion of the Hamiltonian in (19) and (20) beyond the LSWT yields
two anharmonic terms, cubic and quartic, describing three- and four-magnon
interaction, respectively. The cubic anharmonicity comes from the odd part of
the model (20) and carries an important umklapp-like ${\bf Q}$-shift of the
momentum in the one-to-two-magnon coupling, similar to the other models with
the staggered structure of the cubic terms studied in the past; see Refs. [69,
70, 50]. The quartic term is from the even part of the model (19), in which
higher $1/S$-terms of the HP bosonization of spins are kept.

Diagrammatically, these interactions result in a loop expansion, with the
lowest-order diagrams shown in Figs. 5(a)-(c). In a strict $1/S$ sense, their
contributions to the magnon excitation spectrum are of the same order. Three
more diagrams with the same number of loops, corresponding to the anomalous
self-energies, are not shown as they yield corrections of the higher
$1/S$-order [71].

Deferring some essential but technical details concerning three-magnon vertex
symmetrization to Appendix B, the two self-energies in Figs. 5(a) and 5(b) are
the decay and the source ones, respectively,

$\displaystyle\Sigma^{(3)}({\bf k},\omega)=\Sigma^{(d)}({\bf
k},\omega)+\Sigma^{(s)}({\bf k},\omega)\,.$ (27)

While both come from the same cubic anharmonicities, it is the decay diagram
that is relevant to the description of some of the most dramatic modifications
that may occur in the magnon excitation spectra, such as the anomalous
broadening due to quasiparticle breakdown [27, 49, 46], strong renormalization
due to avoided decays [72, 73], and threshold singularities [74, 46]. These
effects occur when the single-particle and two-particle spectra overlap, with
the lower dimensions of the spin system [75], symmetry of the spin model [49],
and favorable kinematics [48, 46, 76] all playing a significant role in the
resultant magnitude of these effects.

All these phenomena manifest themselves quite spectacularly in the CoNb2O6
excitation spectrum in the transverse-field-induced polarized phase [42, 43],
owing to the 1D nature of the zigzag chains, low spin-symmetry leading to a
direct one-to-two-magnon coupling (20), and a favorable overlap with the ${\bf
Q}$-shifted two-magnon continuum, also allowing for the field-variation of it,
the features thoroughly discussed in Ref. [43].

Therefore, analytical insights by the $1/S$ nonlinear SWT (NLSWT) into the
magnon interactions can be expected to shed further light on the important
aspects of the decays, level repulsion, and singularities in the excitation
spectra of CoNb2O6, related zigzag chain materials, and other anisotropic-
exchange magnets. However, this expectation is undermined by the unphysical
divergences in NLSWT at the critical field, characteristic to anisotropic
models [50, 53].

The problem can be seen in the strict $1/S$-expansion for the magnon spectrum,
in which corrections to the LSWT energy (25) are given by the on-shell
($\omega\\!=\\!\varepsilon_{\bf k}$) self-energies from Figs. 5(a)-(c)

$\displaystyle\tilde{\varepsilon}_{\bf k}=\varepsilon_{\bf
k}+\delta\varepsilon_{\bf k}^{(3)}+\delta\varepsilon_{\bf
k}^{(4)},\quad\Gamma_{\bf k}=-\text{Im}\big{[}\Sigma^{(3)}({\bf
k},\varepsilon_{\bf k})\big{]},\quad$ $\displaystyle\delta\varepsilon_{\bf
k}^{(3)}=\text{Re}\big{[}\Sigma^{(3)}({\bf k},\varepsilon_{\bf
k})\big{]},\quad\delta\varepsilon_{\bf k}^{(4)}=\Sigma^{HF}({\bf k}),$ (28)

where $\tilde{\varepsilon}_{\bf k}$ is the renormalized spectrum, $\Gamma_{\bf
k}$ is the decay-induced broadening, $\Sigma^{(3)}({\bf k},\omega)$ from (27)
is discussed above, and the $\omega$-independent Hartree-Fock self-energy
$\Sigma^{HF}({\bf k})$ is shown Fig. 5(c), see Appendix B.3.

Our Figure 6 shows the NLSWT result of such a $1/S$-renormalization of the
magnon spectrum (28), calculated using the best-fit model for CoNb2O6 from
(12) and (15) and for the field just above the classical value of the critical
one in (17), $H\\!=\\!1.01H_{c}$, all for $S\\!=\\!1/2$. An artificial
broadening of $10^{-3}$ meV was used in calculating $\Sigma^{(3)}({\bf
k},\varepsilon_{\bf k})$ (27). Also shown are the LSWT single-magnon branch
(25) together with the ${\bf Q}$-shifted two-magnon LSWT continuum, black
dashed line and the shaded area, respectively.

As is expected, significant singular modifications of the spectrum close to
the decay threshold boundaries, which correspond to the crossing of the
single-magnon branch with the edges of the two-magnon continuum, are present
in the NLSWT spectrum in Fig. 6. The decay-induced scattering rate
$\Gamma_{\bf k}$, divergent at the same thresholds ${\bf k}^{*}$ and ${\bf
G}-{\bf k}^{*}$, with ${\bf G}\\!=\\!2{\bf Q}$, is also shown. This is in
accord with the similarly stark modifications of the magnon spectra in a
variety of other models [71, 70]. Not only do they demonstrate the anomalous
broadening and strong repulsion of the single-magnon spectrum from the two-
magnon continuum within the limited capacity of the naïve perturbation theory,
but they also signify a breakdown of the $1/S$-expansion in the vicinity of
the decay thresholds and call for a more consistent treatment of these
effects, going beyond the $1/S$-approximation to regularize the associated
divergences [71, 46, 49, 48]. We offer further analysis of these physical
threshold singularities in Sec. IV.

Figure 6: Magnon spectrum in the LSWT (25) and NLSWT (28)
$1/S$-approximations, two magnon continuum (shaded area), and $1/S$ decay rate
$\Gamma_{\bf k}$ for $H\\!=\\!1.01H_{c}$, best-fit parameters for CoNb2O6, and
$S\\!=\\!1/2$. Arrows indicate the unphysical divergences and physical
threshold singularities, see the text.

However, this discussion of the physical aspects of the interacting magnon
spectrum is completely undermined, as the results in Fig. 6 are dominated
instead by the nearly divergent behavior of the spectrum near the
$\bm{\Gamma}$ point, which is far away from the decay thresholds. Therefore,
it should not be affected by the anharmonicities and should correspond to a
minimum of the magnon mode in the ferromagnetically-dominated model.

This behavior is due to the $1/S$-expansion, which can also be seen as an
expansion in $1/\varepsilon_{\bf k}$. Because of the closing of the excitation
gap in (26), the $1/S$-corrections to the magnon energy in (28) diverge as
$1/|{\bf k}|$ at $H\\!\rightarrow\\!H_{c}$.

While clearly unphysical, this failure of the NLSWT in the proximity of the
field-induced transition with vanishing excitation gap is not unexpected, as
it is characteristic of the models with the lower spin symmetry, such as
anisotropic-exchange ones [50, 53].

In order to analyze the technical anatomy of this failure, it is instructive
to consider the Hartree-Fock self-energy in Fig. 5(c). To derive it from the
quartic terms in (19), one decouples the four-boson combinations from the
$1/S$-expansion down to the two-boson ones using the real-space HF averages
$\\{\langle a^{\dagger}_{i}a_{i}\rangle,\langle
a^{\dagger}_{i}a_{j}\rangle,\langle
a^{\dagger}_{i}a_{j}^{\dagger}\rangle,\dots\\}$, which can be
straightforwardly evaluated from the Bogolyubov parameters of the LSWT; see
Appendix B for the explicit expressions and technical steps [70, 71].

As a result, the quartic terms are reduced to the LSWT form of Eq. (22), only
with the $1/S$-corrections $\delta A_{\bf k}$ and $\delta B_{\bf k}$ instead
of the $A_{\bf k}$ and $B_{\bf k}$ terms in (23). Then, the $1/S$-expansion
yields

$\delta\varepsilon_{\bf k}^{(4)}=\Sigma^{HF}({\bf k})=\frac{A_{\bf k}\delta
A_{\bf k}-B_{\bf k}\delta B_{\bf k}}{\varepsilon_{\bf k}},$ (29)

which is explicitly divergent in $1/\varepsilon_{\bf k}$ at $H_{c}$ due to the
gapless mode at the $\bm{\Gamma}$ point. One should note that in the highly-
symmetric spin-isotropic models such an expansion is benign, because $\delta
A_{\bf k}$ and $\delta B_{\bf k}$ follow the same ${\bf k}$-dependence as the
LSWT $A_{\bf k}$ and $B_{\bf k}$ terms, canceling the divergence for the
vanishing $\varepsilon_{\bf k}$ [56, 77, 78].

It is important to observe that the strong divergence in Eq. (29) originates
from the strict use of the $1/S$-approximation that can be straightforwardly
avoided by replacing $A_{\bf k}\\!\rightarrow\\!A_{\bf k}+\delta A_{\bf k}$
and $B_{\bf k}\\!\rightarrow\\!B_{\bf k}+\delta B_{\bf k}$ in the renormalized
spectrum, an approach used in a variety of models [50, 79]. In our case, the
solution is more subtle, first because of the cubic terms, but also because of
the 1D character of the problem, which leads to the logarithmically divergent
real-space HF averages for the gapless spectrum. Nevertheless, such an
approach hints at the self-consistent regularization scheme, which does not
only remove the singularity at $H\\!\rightarrow\\!H_{c}$, but also allows us
to access the field range that is inaccessible to the standard spin-wave
theory. This method is discussed next.

###  III.4 Self-consistent Hartree-Fock method

For the CoNb2O6 model discussed in this work, there is a clear hierarchy of
the exchange terms, with the dominant Ising exchange $J_{z_{0}z_{0}}$; see Eq.
(12). Moreover, since the staggered $J_{y_{0}z_{0}}$ term enters only via the
higher-order anharmonic coupling, one can expect that its contribution to the
magnon spectrum away from the threshold singularities is perturbatively small,
$\Sigma^{(3)}({\bf k},\omega)\\!\sim\\!{\cal
O}(J_{y_{0}z_{0}}^{2}/J_{z_{0}z_{0}})$, while the role of the problematic
correction from the quartic terms is $\Sigma^{HF}({\bf k})\\!\sim\\!{\cal
O}(J_{z_{0}z_{0}})$. This suggests the following two-step regularization
procedure.

At the first stage, we neglect the contribution of the cubic terms and perform
a self-consistent calculation of the renormalized eigenvalues
$\bar{\varepsilon}_{\bf k}$ and eigenstates $\bar{u}_{\bf k}$ and
$\bar{v}_{\bf k}$ of the SWT using an iterative procedure with the quartic-
term contribution, referred to as the self-consistent Hartree-Fock (SCHF)
method. It goes beyond the standard SWT by combining different orders in $1/S$
[54, 55, 56], as is depicted in Fig. 5(d), which emphasizes the self-
consistency in the inner line of the HF self-energy. The self-consistency loop
is depicted below

$\displaystyle\begin{array}[]{ccc}\\{\mbox{HFs}\\}&\Longrightarrow&\\{\delta\bar{A}_{\bf
k},\delta\bar{B}_{\bf k}\\}\\\ \Uparrow&&\Downarrow\\\
\\{\bar{\varepsilon}_{\bf k},\bar{u}_{\bf k},\bar{v}_{\bf
k}\\}&\Longleftarrow&\\{\bar{A}_{\bf k},\bar{B}_{\bf k}\\}\end{array}$ (33)

The set of the real-space HF averages, denoted as $\\{\mbox{HFs}\\}$, is used
to obtain the quartic-term contributions to the harmonic theory,
$\delta\bar{A}_{\bf k}$ and $\delta\bar{B}_{\bf k}$, as is described in
Appendix B.3. They result in the modified, but an LSWT-like eigenvalue problem
of the same form as in Eq. (22) with $\bar{A}_{\bf k}\\!=\\!A_{\bf
k}+\delta\bar{A}_{\bf k}$ and $\bar{B}_{\bf k}\\!=\\!B_{\bf
k}+\delta\bar{B}_{\bf k}$, which, in turn, leads to the new set of the
energies $\bar{\varepsilon}_{\bf k}$ and Bogolyubov parameters $\bar{u}_{\bf
k}$ and $\bar{v}_{\bf k}$, with the latter used as an input for the HF
averages. The cycle is continued until numerical convergence in the HF
averages is reached.

With the additional important technical details and step-by-step
implementation described in Appendix B.4, the following point should be made.
Since the described approach regularizes the energy gap in the magnon
spectrum, the gap remains finite at the nominal LSWT critical field $H_{c}$.
Because of that, the SCHF method also allows us to extend our study to the
field values below $H_{c}$, the feat which is unattainable by the standard
$1/S$ SWT approaches. For that, we perform the SCHF calculations at a field
$H\\!>\\!H_{c}$, and then use the outcome for the converged HF averages at a
higher field as an input for the next SCHF calculation for a continuously
decreasing field. The stability and consistency of this procedure is verified
by varying the initial field, the step in the field decrease, and other
iterative parameters; see Appendix B.4.

Finally, the magnon energy spectrum is obtained by reinstating the cubic terms
(27). Importantly, all Bogolyubov coefficients that enter decay and source
vertices as well as the magnon energies in the loops of the diagrams in Figs.
5(a) and 5(b) (see Appendix B.3) are replaced with their regularized values
obtained within the SCHF method described above. Then, the regularized on-
shell ($\omega\\!=\\!\bar{\varepsilon}_{\bf k}$) magnon energy is given by

$\displaystyle\tilde{\varepsilon}_{\bf k}=\bar{\varepsilon}_{\bf
k}+\delta\bar{\varepsilon}_{\bf k}^{(3)},$ (34)
$\displaystyle\delta\bar{\varepsilon}_{\bf
k}^{(3)}=\text{Re}\big{[}\Sigma^{(3)}({\bf k},\bar{\varepsilon}_{\bf
k})\big{]},\quad\bar{\Gamma}_{\bf k}=-\text{Im}\big{[}\Sigma^{(3)}({\bf
k},\bar{\varepsilon}_{\bf k})\big{]}.$

The success of the described regularization procedure in removing the
unphysical singularity at the critical field and in describing the physical
spectrum of CoNb2O6 is demonstrated in the next Section. We also note that in
the results for the dynamical structure factor discussed below, the full
$\omega$-dependence of the cubic self-energies in (27) is used. In the
following, we refer to the method outlined here as the SCHF$+\Sigma^{(3)}$.

##  IV Results

In this Section we present the outcome of the self-consistent method advocated
above, demonstrating its power in regularizing the unphysical divergences and
ability to extend the theory beyond the restrictive classical boundaries. We
also compare our results for the dynamical structure factor in the field-
induced paramagnetic phase with the inelastic neutron scattering data of Refs.
[42, 43].

###  IV.1 Regularization of the unphysical divergences

Figure 7: Same as in Fig. 6. LSWT (25), NLSWT (28), and SCHF$+\Sigma^{(3)}$
(34) magnon spectra for $H=1.01H_{c}$.

The success of our method (34) in regularizing the unphysical divergences
discussed in Sec. III.3 is demonstrated in Fig. 7, where the magnon energy
spectrum by the SCHF$+\Sigma^{(3)}$ for $H\\!=\\!1.01H_{c}$ and the best-fit
model of CoNb2O6 is shown together with the LSWT (25) and NLSWT (28) results
from Fig. 6. In the regularized spectrum in Fig. 7, the offending divergent
behavior of the NLSWT near the ${\bm{\Gamma}}$ point is gone, and one is able
to focus on the physical effects of magnon interaction in the decay-related
phenomena.

The second achievement is the following. The 1D critical field calculated by
DMRG for the CoNb2O6 model, $B_{c}^{\text{1D}}\\!\approx\\!4.52(1)\
\mathrm{T}$, is close to the experimentally estimated one,
$B_{c}^{\text{1D}}\\!=\\!5\ \mathrm{T}$ [40], both much smaller than the
classical critical field, $B_{c}^{\text{cl}}\\!\approx\\!8.8\ \mathrm{T}$,
obtained for the best-fit parameters in Eq. (17). For
$B\\!<\\!B_{c}^{\text{cl}}$, the paramagnetic phase is not a minimum of the
classical energy and cannot be studied by means of the $1/S$-expansion,
because the LSWT Hamiltonian (22) is not positive-definite, with its spectrum
(25) becoming imaginary in some regions of ${\bf k}$.

Our Fig. 8 shows the magnon excitation gap $\Delta_{0}$ at the ${\bm{\Gamma}}$
point as a function of the field for the best-fit model of CoNb2O6 obtained by
different methods. The vanishing of this gap corresponds to a phase transition
from the paramagnetic to the ordered phase at $T\\!=\\!0$. The red horizontal
line on top of the figure and the gray shaded area emphasize the difference
between the classical and experimental results for it. According to the LSWT,
the gap vanishes at the classical critical field (17), and it diverges in the
NLSWT $1/S$-approximation. The SCHF$+\Sigma^{(3)}$ method regularizes this
divergence at $B_{c}^{\text{cl}}$ and allows us to extend the study of the
magnon spectrum into the field region below the classical boundary to the
paramagnetic phase. These results also shows an excellent agreement with the
experimental data for the gap at 7, 8, and 9 $\mathrm{T}$ [80], highlighting
the quantitative accuracy of our approach.

Last but not the least, Fig. 8 shows the results of the DMRG simulations for
the gap in the same model, which agree closely with both experimental data and
results of the self-consistent theory, except for the close proximity of the
critical field. With the details of the DMRG calculations deferred to Appendix
C, it should be noted that the DMRG critical field for the single-chain 1D
model is $B_{c}^{\text{1D}}\\!=\\!4.52(1)\ \mathrm{T}$, also below the
experimental value. This is because the 3D interchain terms play important
role near the transition; see Refs. [44, 45].

Deviation of the SCHF$+\Sigma^{(3)}$ from the DMRG results is of a different,
but related nature. While remarkably successful otherwise, the self-consistent
method is not entirely consistent for the gap approaching zero. The SCHF
method naturally prevents the gap from closing, because the HF averages would
diverge logarithmically for the gapless spectrum due to the 1D character of
the model. The finite critical field in the SCHF$+\Sigma^{(3)}$ method is due
to the non-self-consistent perturbative treatment of the cubic term
$\Sigma^{(3)}$, and the value of the critical field it yields at
$\approx\\!3.4\ \mathrm{T}$ is not physically meaningful.

Notwithstanding these minor limitations and concerns, one should not lose the
sight offered by Fig. 8, which demonstrates the ability of our approach to
provide quantitatively meaningful description of the magnetic excitations for
a wide field range in the paramagnetic phase of CoNb2O6. This is despite
strong quantum fluctuations, anisotropic exchanges, and low dimensionality of
the problem, the factors that make the standard SWT fail. Not only does the
SCHF$+\Sigma^{(3)}$ method regularize the unphysical divergences, but it also
preserves the physical features of the threshold singularities, which are
discussed next.

Figure 8: Magnon gap $\Delta_{0}$ vs $B$ by LSWT (25), NLSWT (28),
SCHF$+\Sigma^{(3)}$ (34), and DMRG for the best-fit model of CoNb2O6, compared
to experimental gaps (solid points) [80]. The DMRG and classical critical
fields and the region inaccessible by standard SWT (gray shaded area) are
highlighted, see text.

###  IV.2 Decay threshold singularities

While the main results and comparison with the experimental data for the
dynamical structure factor will be discussed in the next section, we would
like to briefly recall the origin and the nature of the decay-related
phenomena in the magnon spectra; see also Refs. [69, 46, 71].

In Fig. 9 we show the magnon spectrum $\bar{\varepsilon}_{\bf k}$ obtained by
SCHF method discussed in Sec. III.4, together with $\tilde{\varepsilon}_{\bf
k}$ from Eq. (34) that includes contribution of the on-shell cubic self-energy
$\Sigma^{(3)}({\bf k},\bar{\varepsilon}_{\bf k})$, for the best-fit model of
CoNb2O6 and $H\\!=\\!0.8H_{c}$ ($B\\!\approx\\!7$ T), well below the classical
critical field $H_{c}$. The two-magnon density of states (DoS), $D^{(2)}({\bf
k},\omega)\\!=\\!\frac{1}{N}\sum_{\bf q}\delta(\omega-\bar{\varepsilon}_{{\bf
q}}-\bar{\varepsilon}_{{\bf k}-{\bf q}+{\bf Q}})$, for the SCHF energies
$\bar{\varepsilon}_{\bf q}$, is shown as an intensity plot.

Figure 9: Magnon spectrum by the SCHF (34) with and without the cubic self-
energy $\Sigma^{(3)}({\bf k},\bar{\varepsilon}_{\bf k})$ for the best-fit
model of CoNb2O6 and $B\\!=\\!7$ T ($H\\!=\\!0.8H_{c}$). The intensity plot is
the two-magnon DoS calculated using the SCHF energies $\bar{\varepsilon}_{\bf
q}$.

As one can see in Fig. 9, the contribution of the cubic term
$\Sigma^{(3)}({\bf k},\bar{\varepsilon}_{\bf k})$ to the magnon spectrum away
from the crossings with the two-magnon continuum is indeed small, as is
anticipated in the discussion of the SCHF$+\Sigma^{(3)}$ approach in Sec.
III.4. In fact, the cubic term vanishes entirely at ${\bf k}\\!=\\!{\bf Q}$,
owing to the staggered structure of the corresponding spin-exchange terms
(20), which, in turn, is translated into the antisymmetric structure of the
cubic vertices, see Appendix B.3.

The divergent behavior exhibited by the on-shell SCHF$+\Sigma^{(3)}$ spectrum
near the crossing with the two-magnon continuum at ${\bf k}^{*}$ and
equivalent points, referred to as the decay threshold boundaries [46], is due
to a resonance-like coupling of the single-magnon branch with the two-magnon
continuum, which is provided by the cubic terms. Since the lowest two-magnon
energy must necessarily correspond to a minimum of $E_{{\bf k},{\bf
q}}\\!=\\!\bar{\varepsilon}_{\bf q}+\bar{\varepsilon}_{{\bf k}-{\bf q}+{\bf
Q}}$ at any given ${\bf k}$, it follows that the corresponding two-magnon DoS
must be singular at that minimum in 1D, as one can observe in Fig. 9; see also
Appendix C.2. It also follows that one can expand the denominator of the decay
part of the on-shell self-energy, $\bar{\varepsilon}_{\bf k}-E_{{\bf k},{\bf
q}}$, near such a minimum in the proximity of the threshold ${\bf k}^{*}$ for
small $\Delta{\bf k}\\!=\\!{\bf k}-{\bf k}^{*}$. Because the decay vertex has
no symmetry constraints at a generic ${\bf k}^{*}$ and must generally be
finite, one can obtain the asymptotic behavior for the real and imaginary
parts of the on-shell self-energy on the two sides of the threshold

$\displaystyle\text{Re}\big{[}\Sigma^{(3)}({\bf k},\bar{\varepsilon}_{\bf
k})\big{]}\propto\begin{cases}-1/\sqrt{|\Delta{\bf k}|},&\text{for}\
\Delta{\bf k}<0,\\\ \Lambda+\gamma\Delta{\bf k},&\text{for}\ \Delta{\bf
k}>0,\end{cases}$ $\displaystyle\phantom{\text{Re}\big{[}\Sigma^{(3)}({\bf
k},\bar{\varepsilon}\,}\bar{\Gamma}_{\bf k}\propto\begin{cases}0,&\text{for}\
\Delta{\bf k}<0,\\\ 1/\sqrt{\Delta{\bf k}},&\text{for}\ \Delta{\bf
k}>0,\end{cases}$ (35)

where $\Lambda\\!>\\!0$ is the cut-off parameter and $\gamma$ is a constant.
The inverse square-root singularities in (35) are from the 1D Van Hove
singularity at the edge of the two-magnon continuum that gets imprinted on the
single-magnon spectrum via the anharmonic coupling. These asymptotic results
explain the behavior observed in Fig. 9 and should be contrasted with a
typically weaker singularities in the higher dimensions and for the more
symmetric models [46, 71].

Note that, given the relative simplicity of the magnon dispersion in the
paramagnetic phase of CoNb2O6, the two-magnon energy, $\bar{\varepsilon}_{\bf
q}+\bar{\varepsilon}_{{\bf k}-{\bf q}+{\bf Q}}$, can be well-approximated as
an energy of two particles with the nearest-neighbor hopping,
$\bar{\varepsilon}_{\bf q}\\!\approx\\!{\sf E_{0}}+{\sf J}_{1}\gamma_{\bf
q}^{(1)}$, total momentum ${\bf k}$, and the ${\bf Q}$ shift,
straightforwardly yielding the bow-tie form of the continuum with zero width
at the ${\bm{\Gamma}}$ point, see Fig. 9. In that case, the minimum of the
two-magnon energy corresponds to the energy of two magnons with equivalent
momenta, $E^{\rm min}_{{\bf k},{\bf q}^{*}}\\!=\\!2\bar{\varepsilon}_{{\bf
q}^{*}}$ with ${\bf q}^{*}\\!=\\!({\bf k}-{\bf Q})/2$, for any ${\bf k}$. This
latter condition also holds for most ${\bf k}$ for the true form of
$\bar{\varepsilon}_{\bf q}$, see Appendix C.2.

However, because of the further-neighbor exchanges and a relativistic form of
the magnon dispersion, there is more structure in the two-magnon continuum in
the vicinity of the ${\bm{\Gamma}}$ point, providing the bow-tie region with a
finite width and a richer set of the Van Hove singularities visible in Fig. 9.
Since they are far away from the physical decay thresholds, we refrain from
discussing them here and delegate a more detailed analysis of the two-magnon
kinematics to Appendix C.2.

As is discussed above, the perturbative consideration of the decay diagram
within the on-shell approach offered by Fig. 9, Eq. (34), and Eq. (35)
signifies a breakdown of the perturbation theory in the vicinity of the decay
thresholds and calls for a more consistent treatment of these effects to
regularize the associated divergences. Qualitatively, upon a self-consistent
treatment of the higher-order contributions, which is typically difficult to
implement, these singularities are expected to lead to the anomalous
broadening, renormalization of the magnon spectrum, and the so-called
termination points [74, 71].

One of the less-sophisticated regularizations, which, nevertheless, avoids the
divergences of the on-shell approach, is based on a straightforward use of the
explicit $\omega$-dependence in the cubic self-energy in (27). It corresponds
to the one-loop approximation for the single-magnon Green’s function $G({\bf
k},\omega)$ and its spectral function $A({\bf k},\omega)$, which is able to
yield the quantitatively faithful description of the quasiparticle-like and
incoherent parts of the single-magnon spectrum [71, 75, 72]. Since the
spectral function is directly related to the dynamical structure factor
$S({\bf k},\omega)$, measured in the inelastic neutron-scattering experiments,
we will use this approach as is discussed in the next Section.

###  IV.3 Dynamical structure factor

The general form of the dynamical structure factor (DSF) for the neutron
scattering is [44]

${\cal S}({\bf
k},\omega)=\sum_{\alpha,\beta}g_{\alpha}g_{\beta}\left(\delta_{\alpha,\beta}-\frac{k_{\alpha}k_{\beta}}{k^{2}}\right){\cal
S}^{\alpha\beta}({\bf k},\omega),$ (36)

with the momentum and energy transfer ${\bf k}$ and $\omega$, axes of the
reference frame $\alpha$ and $\beta$, g-tensor components $g_{\alpha}$, and
the spin-spin dynamical correlation function

$\mathcal{S}^{\alpha\beta}({\bf
k},\omega)=\frac{1}{\pi}\,\text{Im}\int_{-\infty}^{\infty}dt\,e^{i\omega
t}\,i\big{\langle}\mathcal{T}S_{\bf k}^{\alpha}(t)S_{-{\bf
k}}^{\beta}(0)\big{\rangle}\,.$ (37)

For the field-induced fluctuating paramagnetic state of CoNb2O6, with the
choice of the local $\\{x,y,z\\}$ or laboratory $\\{x_{0},y_{0},z_{0}\\}$ axes
in Fig. 2(c), only diagonal components of (37) are considered [44]. Two of
them are in the $ac$ plane normal to the field, one along the Ising axis,
${\cal S}^{xx}$ $({\cal S}^{z_{0}z_{0}})$, and one perpendicular to it, ${\cal
S}^{yy}$ $({\cal S}^{x_{0}x_{0}})$ [41]. One more component is along the
field, ${\cal S}^{zz}$ $({\cal S}^{y_{0}y_{0}})$.

In the studies of the CoNb2O6 spectrum in the polarized phase that are
discussed in Refs. [41, 43, 44], the momentum transfer is not aligned
exclusively along the chain direction $c$, see Fig. 2(a), but has other
components. In our consideration, which is focused on the single spin-chain
model, the additional momentum component along the $b$ axis is important
because it is able to detect the zigzag structure of the chain.

Assuming the momentum transfer in the $bc$ plane and keeping only diagonal
components in the DSF, the general expression in Eq. (36) is simplified to

$\displaystyle{\cal S}({\bf k},\omega)$
$\displaystyle=g_{z_{0}}^{2}\left(1-\tilde{\lambda}\,\cos^{2}\gamma\right){\cal
S}^{xx}({\bf k},\omega)$ (38)
$\displaystyle+g_{x_{0}}^{2}\left(1-\tilde{\lambda}\,\sin^{2}\gamma\right){\cal
S}^{yy}({\bf k},\omega)+g_{y_{0}}^{2}\tilde{\lambda}\,{\cal S}^{zz}({\bf
k},\omega),$

where we use the shorthand notation
$\tilde{\lambda}\\!=\\!k_{c}^{2}/(k_{c}^{2}+k_{b}^{2})$ and the angle $\gamma$
is between the Ising and $c$ axes, as before, with the transfer momentum ${\bf
k}$ in the local frame

${\bf k}=k_{b}\hat{\sf b}+k_{c}\hat{\sf c}=k_{b}\hat{z}+k_{c}(\cos\gamma\
\hat{x}-\sin\gamma\ \hat{y}),$ (39)

see Fig. 2. As is discussed in Ref. [43], the non-zero $b$-component of the
momentum in (39) is responsible for the secondary ${\bf Q}$-shifted signal in
the structure factor from the doubling of the unit cell [41, 43, 44]. The
general form of the diagonal components of the dynamical correlation function
is given by

$\displaystyle{\cal S}^{\alpha\alpha}({\bf k},\omega)=$
$\displaystyle\cos^{2}(k_{b}{\sf b})\,\widetilde{{\cal S}}^{\alpha\alpha}({\bf
k},\omega)$ $\displaystyle+$ $\displaystyle\sin^{2}(k_{b}{\sf
b})\,\widetilde{{\cal S}}^{\alpha\alpha}({\bf k}+{\bf Q},\omega),$ (40)

where ${\sf b}$ is the width of the zigzag chain in the $b$ direction, see
Fig. 2(a) and Ref. [41], and $\widetilde{{\cal S}}^{\alpha\alpha}({\bf
k},\omega)$ is the correlation function that depends only on the momentum in
the $c$ direction, ${\bf k}\\!=\\!k_{c}\hat{\textbf{{\sf c}}}$. We note that
in Eq. (40) the main signal is associated with the first term and the
secondary, “shadow” signal, with the ${\bf Q}$-shifted one.

The $\widetilde{{\cal S}}^{xx}$ and $\widetilde{{\cal S}}^{yy}$ components of
the structure factor are the transverse ones and can be straightforwardly
related to the single-magnon spectral function [41, 72] as

$\widetilde{{\cal S}}^{\alpha\alpha}({\bf
k},\omega)=\mathcal{F}^{\alpha\alpha}({\bf k})A({\bf k},\omega),$ (41)

where the kinematic formfactors

$\mathcal{F}^{xx}_{\bf k}=\frac{S}{2}\,(\bar{u}_{\bf k}+\bar{v}_{\bf k})^{2},\
\ \ \ \mathcal{F}^{yy}_{\bf k}=\frac{S}{2}\,(\bar{u}_{\bf k}-\bar{v}_{\bf
k})^{2},$ (42)

produce the ${\bf k}$-dependent modulation of the single-magnon spectral peaks
throughout the Brillouin zone.

The DSF can also be expected to exhibit significant decay-related features,
such as incoherent parts of the single-magnon spectrum and strong
renormalizations, due to the cubic self-energy (27) in the single-magnon
spectral function $A({\bf k},\omega)\\!=\\!-\frac{1}{\pi}{\rm Im}[G({\bf
k},\omega)]$, where

$\displaystyle G({\bf k},\omega)=\frac{1}{\omega-\bar{\varepsilon}_{\bf
k}-\Sigma^{(3)}({\bf k},\omega)+i0^{+}}\,,$ (43)

is the Green’s function in the SCHF$+\Sigma^{(3)}$ approach.

The DSF component along the field, $\widetilde{{\cal S}}^{zz}$, corresponds to
the longitudinal fluctuations, which account for the direct two-magnon
continuum contribution to it,

$\displaystyle\widetilde{{\cal S}}^{zz}({\bf k},\omega)$
$\displaystyle=\sum_{\bf q}\mathcal{F}^{zz}_{{\bf q},{\bf
k}}\,\delta(\omega-\bar{\varepsilon}_{\bf q}-\bar{\varepsilon}_{{\bf k}-{\bf
q}}),$ $\displaystyle\mathcal{F}^{zz}_{{\bf q},{\bf k}}$
$\displaystyle=\frac{1}{2}(\bar{u}_{\bf q}\bar{v}_{{\bf k}-{\bf
q}}+\bar{v}_{\bf q}\bar{u}_{{\bf k}-{\bf q}})^{2}.$ (44)

Note that in contrast to the two-magnon continuum in the anharmonic coupling,
this continuum is not umklapp-shifted by the momentum ${\bf Q}$.

Figure 10: Intensity plots of the dynamical structure factor ${\cal
S}(\textbf{k},\omega)$ in the paramagnetic phase of CoNb2O6 for the field
$7\mathrm{T}$ [(a) and (c)] and $9\mathrm{T}$ [(b) and (d)], best-fit
parameters [(a) and (b)] and adjusted $\lambda_{s}$ [(c) and (d)].
Experimental data were adapted from Ref. [43]. Artificial broadening of
$10^{-3}\ \mathrm{meV}$ for the self-energy (27) and $5\\!\times\\!10^{-2}\
\mathrm{meV}$ for the Green’s functions (43) and (44) were used. Dashed lines
show the bottom of the two-magnon continuum and the arrows indicate shadow
modes.

For comparison with experiments in Ref. [42], in Fig. 10 we illustrate our
calculations, where we combined all contributions to the structure factor as
given by Eq. (38), and used a momentum transfer in the $bc$ plane with
$k_{b}/k_{c}\\!=\\!0.3$ to have a finite contribution from the shadow mode. We
also used artificial broadenings of $10^{-3}\ \mathrm{meV}$ in the self-energy
$\Sigma^{(3)}({\bf k},\omega)$ (27) and $5\\!\times\\!10^{-2}$ meV in the
Green’s functions (43) and (44).

Figure 10 displays our main results. It shows the compilations of the DSF
intensity maps by the inelastic neutron scattering, adapted from Ref. [43],
with the DSF intensities obtained by the theoretical SCHF$+\Sigma^{(3)}$
approach described above. Each plot consists of the experimental data in the
left panel, theoretical results in the middle panel, and the two of them
overlaid in the right panel, where we have exploited the symmetry of the
results about the $\Gamma$ point. The comparison is presented for the field
$7~{}\mathrm{T}$ in Figs. 10(a) and 10(c) and for $9~{}\mathrm{T}$ in Figs.
10(b) and 10(d), respectively. The upper row, Figs. 10(a) and 10(b), shows
theoretical results for the best-fit model of CoNb2O6, see Sec. II.5, while
the lower row, Figs. 10(c) and 10(d), has one parameter from that set
modified. The same comparison with the experimental data for $8~{}\mathrm{T}$
is given in Appendix C.3.

One can see that the theoretical results for the best-fit model (upper row of
Fig. 10) already yield if not an ideal, but a close quantitative agreement
with the experimental data on the gap and magnon bandwith with no adjustment
to the parameters. We note that the mismatch in the energies at the ${\bf Q}$
point might in part be related to the effect of the weak 3D interchain
coupling, which affect the experimental data, but are not included in our
model. The lower row of Fig. 10 shows that an even closer agreement can be
reached by a modest change of a single parameter in the model,
$\lambda_{s}\\!=\\!(J_{x_{0}x_{0}}\\!+J_{y_{0}y_{0}})/2J_{z_{0}z_{0}}$, used
in the parametrization of Ref. [44], see also Sec. II.5, to which the maximum
of the magnon band is most sensitive. Changing it from the best-fit value of
$\lambda_{s}\\!=\\!0.251$ to $\lambda_{s}\\!=\\!0.2$ improves the agreement
with our theory. This is not to challenge the comprehensive multi-dimensional
best-fit strategy of Ref. [44], but to highlight, once again, that a self-
consistent approach can turn a theory plagued with unphysical divergences into
a reliable, nearly quantitative tool.

Turning to the other features of the theoretical results shown in Fig. 10, it
is clear that the off-shell $\omega$-dependent cubic self-energy successfully
regularizes the threshold singularities discussed in Sec. IV.2, and, indeed,
provides a quantitatively faithful description of the quasiparticle-like and
incoherent parts of the single-magnon spectrum. On the inner side of the two-
magnon continuum, the magnon spectral lines acquire a substantial broadening
in a close agreement with the experimental data, see also the plot for
$8~{}\mathrm{T}$ in Appendix C.3.

On the outer side of the continuum, a direct intersect of the magnon mode with
the continuum is avoided via a strong renormalization of the magnon energy,
creating a gap-like splitting and a characteristic loss of the spectral weight
of the magnon line. Although the quantitative agreement with the experimental
data on the size of the gap-like feature and its evolution with the field is
rather spectacular, the theoretical results contain more details, with the
remnant of the magnon mode following the bottom of the two-magnon continuum
for an extended range of the momenta. While a recent proposal suggests that in
1D such an edge-mode should survive for all the momenta [75], this conclusion
is an artifact of the one-loop approximation for the magnon self-energy, which
is also employed in our study. In reality, it is expected that the magnon mode
should meet the continuum at the so-called termination point [46, 74, 71], the
result that requires a self-consistent treatment of the higher-order diagrams
in the theory, which is not attempted here.

Because of the constraint provided by the quantitative accord of the
experiment and theory on the size of the gap-like splitting in the spectrum,
an additional comment can be made on the potential value of the “residual”
$J_{x_{0}y_{0}}$ term in the exchange matrix (6), discussed in Sec. II.4.2 and
Appendix A.2. This term does not modify the even part of the spin Hamiltonian
(19), but contributes to the cubic coupling from the odd part (20). Because of
the staggered nature of this term, the structure of the decay vertex is not
expected to modify, leading to an enhancement of the decay self-energy
according to $\Sigma^{(3)}({\bf
k},\omega)\\!\propto\\!J_{x_{0}y_{0}}^{2}\\!+J_{y_{0}z_{0}}^{2}$. Since the
best-fit parameters without the $J_{x_{0}y_{0}}$-term provide a close
quantitative description of the decay-related features described above, see
Fig. 10, one can conclude that the $J_{x_{0}y_{0}}$ term must be small
compared to $J_{y_{0}z_{0}}$. This observation supports the approach of Refs.
[43, 44] and our discussion in Sec. II.4.2.

In Fig. 10, in both theoretical and experimental results, one can also observe
the ${\bf Q}$-shifted “shadow” mode in addition to the main contribution from
the single-magnon excitations [41, 43, 44]. As is discussed above, it
originates from the non-zero component of the transfer momentum along the $b$
axis. In the theory results in Fig. 10, one can also observe a continuum-like
contribution from the longitudinal component of the structure factor (44),
with its role being generally minor.

##  V Longitudinal field effects

Up to this point, our study of the CoNb2O6 excitation spectrum in the
paramagnetic phase concerned the transverse direction of the field. Now we
focus on the effects of an additional longitudinal field component. In this
consideration, we use the same spin Hamiltonian
$\hat{\mathcal{H}}_{1}\\!+\\!\hat{\mathcal{H}}_{2}$, Eqs. (1), (7), and (14),
with the transverse-field term
$\hat{\mathcal{H}}_{\perp}\\!=\\!g_{y_{0}}\mu_{B}B_{\perp}\sum_{i}S_{i}^{y_{0}}$,
Eq. (16), now augmented by the longitudinal-field term

$\hat{\mathcal{H}}_{\parallel}=g_{z_{0}}\mu_{B}B_{\parallel}\sum_{i}S_{i}^{z_{0}}.$
(45)

with the field component $B_{\parallel}$ along the Ising $z_{0}$ axis; see
Fig. 2. We are interested in the regime of the weak longitudinal fields,
$B_{\parallel}\\!\ll\\!B_{\perp}$, with excitations remaining spin-flip-like
and magnon description of the their spectrum still adequate.

Although symmetry-wise the classification of the type of the symmetry-breaking
provided by the longitudinal field in (45) in the case of the zigzag model of
CoNb2O6 is more delicate, relating it to the glide-symmetry breaking [43], the
main effect is the same as in the paradigmatic Ising model [81, 82]. The
quantum phase transition of the transverse-field Ising-like model ceases to
exist and turns into a crossover, with a finite excitation gap at the former
transition point.

The second effect of the symmetry-breaking longitudinal field is specific to
the zigzag chain model, and it is intimately tied to the presence of the
staggered bond-dependent $J_{y_{0}z_{0}}$ terms, allowed by the same glide
symmetry. Because of the tilt of the spin-quantization axis induced by the
longitudinal field in the paramagnetic state, the two-site unit cell of the
zigzag structure becomes explicit in the model of spin flips, doubling the
unit cell of the “simple” Ising chain that sufficed until now. The description
of the excitation spectrum in this more general case requires two distinct
branches of spin excitations within the reduced Brillouin zone of the zigzag
chain.

Importantly, these two excitation branches will be split by a band gap. As we
argue below, one can expect strong modifications of the two-magnon DoS as a
result of these changes in the single-magnon spectrum, inducing richer
varieties of the Van Hove singularities that are potentially observable.
Below, we quantify both effects using the LSWT formalism.

###  V.1 The excitation gap and the band gap

In the tilted field with a small longitudinal component away from the
transverse $y_{0}$ axis toward the Ising $z_{0}$ axis, the spin-quantization
axis will tilt by the angle $\theta$ in the $y_{0}z_{0}$ $(zx)$ plane, see
Fig. 2(c), found from the minimization of the classical energy, see Appendix D
for details,

$H_{\perp}\sin\theta-H_{c}\sin\theta\cos\theta-H_{\parallel}\cos\theta=0,$
(46)

where $H_{\perp}\\!=\\!g_{y_{0}}\mu_{B}B_{\perp}$ and
$H_{\parallel}\\!=\\!g_{z_{0}}\mu_{B}B_{\parallel}$ are the transverse and
longitudinal fields, respectively, in the energy units, and the classical
critical field $H_{c}$ is from Eq. (17).

Here, we focus on the case of the transverse field value equal to the critical
field, $H_{\perp}\\!=\\!H_{c}$, and study the dependence of the spectrum gap
and the band gap on the longitudinal field, as both gaps vanish at
$H_{\parallel}\\!=\\!0$. For $H_{\parallel}\\!\ll\\!H_{\perp}\\!=\\!H_{c}$,
Eq. (46) can be solved by expanding in the small canting angle, yielding

$\theta=\left(2H_{\parallel}/H_{c}\right)^{1/3}.$ (47)

Note that this fractional power law is reminiscent of that of the canting
angle due to staggered Dzyaloshinskii-Moriya interaction in the isotropic
square-lattice antiferromagnet near saturation field [83].

The LSWT consideration of the $1/S$-expansion of the model with the
longitudinal field involves a somewhat cumbersome diagonalization of the
$4\times 4$ Hamiltonian for the two bosonic species, deferred to Appendix D,
which gives explicit expressions of the energies of the two magnon branches.
In the presence of the longitudinal field, excitation gap $\Delta_{0}$ at the
${\bf\Gamma}$ point and the band gap $\Delta_{b}$ at the ${\bf Q}/2$ point
open up, as is discussed above; see Fig. 11(a) for a comparison to the
$H_{\parallel}\\!=\\!0$ case.

Figure 11: LSWT results for the best-fit model of CoNb2O6 and
$H_{\perp}\\!=\\!H_{c}$. (a) Magnon spectrum for $H_{\parallel}\\!=\\!0$
(dashed line) and $H_{\parallel}\\!=\\!0.02H_{c}$ (solid line), with gaps
$\Delta_{0}$ and $\Delta_{b}$ identified. (b) Excitation gap $\Delta_{0}$ and
band gap $\Delta_{b}$ vs $(H_{\parallel}/H_{c})^{1/3}$. Lines are the full
LSWT results (Appendix D) and dashed lines are their asymptotics from Eq.
(48).

Using the canting angle (47) for small fields, one can obtain asymptotic
expressions for the gaps,

$\displaystyle\Delta_{0}\approx\alpha_{0}H_{\parallel}^{1/3},\ \
\alpha_{0}\approx\sqrt{\frac{3}{2}}H_{c}\left(\frac{2}{H_{c}}\right)^{1/3},$
$\displaystyle\Delta_{b}\approx\alpha_{b}H_{\parallel}^{1/3},\ \
\alpha_{b}\approx
2S\big{|}J_{y_{0}z_{0}}\big{|}\left(\frac{2}{H_{c}}\right)^{1/3},$ (48)

which follow the same fractional power law vs field, see Appendix D for the
exact proportionality coefficients and Fig. 11(b), which shows a comparison of
the asymptotic results (48) with the full LSWT results for the best-fit model
of CoNb2O6.

We note that, according to the Ising conformal field theory in $1+1$
dimensions [81], the scaling of the spectrum gap with the longitudinal field
is known to obey a different fractional power law with the exponent $8/15$.
Still, the expressions in Eq. (48) highlight an important distinction of the
two gaps. The spectrum gap $\Delta_{0}$ is essentially the same as it would
have been for the “simple” Ising-like spin chain, as it is independent of the
bond-dependent terms. However, the appearance of the band gap $\Delta_{b}$ is
precisely due to the staggered bond-dependent $J_{y_{0}z_{0}}$ terms, rooted
in the zigzag nature of the model.

Quantitatively, because the bond-dependent terms in CoNb2O6 are secondary to
the main Ising term, the excitation gap $\Delta_{0}$ grows faster with the
longitudinal field than the band gap $\Delta_{b}$. Comparison of their
asymptotics in Eq. (48) for the CoNb2O6 model yields

$\frac{\Delta_{b}}{\Delta_{0}}\approx
2S\sqrt{\frac{2}{3}}\,\frac{\big{|}J_{y_{0}z_{0}}\big{|}}{H_{c}}\approx 0.27.$
(49)

As we will see next, this result has a significant impact on the longitudinal
field range for which the overlap of the single-magnon branches with the
additional Van Hove singularities in the two-magnon spectra is possible.

###  V.2 More threshold singularities

One of the important consequences of the magnon band splitting is the explicit
separation of the two-magnon continuum into three continua, corresponding to
different combinations of the single-magnon species

$E_{{\bf k},{\bf q}}^{\\{\mu,\nu\\}}=\varepsilon_{\mu{\bf
q}}+\varepsilon_{\nu{\bf k}-{\bf q}},$ (50)

where $\mu(\nu)\\!=\\!1,2$. This splitting also necessarily creates richer
structure of the Van Hove singularities in the continuum, which can affect the
single-magnon spectrum via the anharmonic coupling. Thus, if allowed by the
two-magnon kinematics, the longitudinal field can potentially lead to more
singularities in the magnon spectra, in addition to the ones discussed in
Secs. IV.2 and IV.3.

In Fig. 12, we show magnon spectrum together with the two-magnon DoS intensity
plot for $H_{\perp}\\!=\\!H_{c}$ and $H_{\parallel}\\!=\\!0.01H_{c}$
($B_{\parallel}\\!\approx\\!0.09\ \mathrm{T}$) for the best-fit model of
CoNb2O6, from which one can appreciate the more intricate structure of the
field-induced Van Hove singularities in the two-magnon continuum.

However, in practice, because excitation gap $\Delta_{0}$ grows faster than
the band gap $\Delta_{b}$ (49), such a trend in the field-induced gaps
provides a rather narrow range of the longitudinal fields for which the
kinematics is favorable of the crossing of the additional Van Hove
singularities by the single-magnon spectrum. Thus, already for
$H_{\parallel}\\!=\\!0.01H_{c}$ shown in Fig. 12, the magnon branch barely
accesses the extra features in the continuum, prohibiting such a crossing for
the larger fields.

In addition to the kinematics, we would also like to remark on the effect of
the small longitudinal field on the structure of the anharmonic cubic term
that is responsible for the one-to-two-magnon coupling. There are two parts in
it in the presence of the $H_{\parallel}$ field, one that largely retains the
same form as in the odd part of the Hamiltonian (20), originating from the
staggered exchanges $J_{y_{0}z_{0}}$, while the other is due to the tilt angle
$\theta\\!\propto\\!H_{\parallel}^{1/3}$ that allows most other exchanges to
contribute to the cubic anharmonicity; see Appendix D for some more detail. We
note that it is the latter term that was previously considered as the main
source of the decay singularities in CoNb2O6 [42]. However, not only it is
subleading in the weak longitudinal-field regime, but it is also not
staggered, resulting in an unfavorable kinematics for the decay-related
processes.

Figure 12: Magnon spectrum for $H_{\perp}\\!=\\!H_{c}$,
$H_{\parallel}\\!=\\!0.01H_{c}$ and the best-fit model of CoNb2O6 in the
repeated Brillouin zone scheme. The intensity plot is the two-magnon DoS for
the continua in (50).

##  VI Conclusions

We conclude by summarizing our results. In this study, we have thoroughly
reanalyzed the symmetry-based approach to formulating anisotropic-exchange
model of the quasi-one-dimensional ferromagnet CoNb2O6. We have proposed a
connection of its model to the broader class of such models, studied for a
wide variety of materials with complex bond-dependent spin-orbit-induced
exchanges. We have also clarified the role of a phenomenological constraint
that has been used to restrict parameter space of CoNb2O6, and have
investigated the magnitude and the effects of the residual terms, which were
neglected in the previous studies, using real-space perturbation theory and
unbiased DMRG approach.

The main result of the present work is the self-consistent study of the
effects of magnon interactions in the excitation spectrum of CoNb2O6 in the
quantum paramagnetic phase. We have proposed and applied a self-consistent
Hartree-Fock regularization of the problematic unphysical divergences in the
$1/S$ spin-wave expansion that is common to various anisotropic-exchange
models. Not only does this method eliminate such unphysical singularities, but
it also preserves the integrity of the threshold phenomena of magnon decay and
spectrum renormalization that are present in both theory and experiment of
CoNb2O6.

Using the microscopic parameters proposed previously, we have employed this
approach to study excitation spectrum in the fluctuating paramagnetic phase of
CoNb2O6. For the dynamical structure factor ${\cal S}({\bf k},\omega)$, we
have demonstrated a close quantitative agreement of our theory with the
neutron-scattering data for both the quasiparticle-like and incoherent parts
of the single-magnon spectrum, also in the field regime that is inaccessible
by the standard spin-wave theory. Moreover, our results for the spectrum gap
are in a close accord with the complimentary DMRG calculations for the same
model parameters.

These results prove the ability of our approach to provide quantitatively
faithful description of the magnetic excitations in the paramagnetic phase of
CoNb2O6, despite strong quantum fluctuations, anisotropic exchanges, and low
dimensionality of the problem, the factors that make the standard SWT-like
approaches fail. Furthermore, it can be expected that our approach should be
able to yield further analytical insights into magnon interactions and decay
phenomena and shed light on the important aspects of the excitation spectra in
the other anisotropic-exchange magnets in the higher dimensions, where it is
free from the remaining minor inconsistencies associated with the 1D nature of
the CoNb2O6 model.

Lastly, we have also discussed the effects of additional longitudinal fields
in the paramagnetic phase of CoNb2O6. We have demonstrated that due to the
zigzag lattice structure and affiliated bond-dependent exchanges in the model,
and due to the symmetry breaking by the longitudinal field, both the
excitation gap and the band gap develop in the magnon spectrum. We have
described how the band splitting leads to the additional anomalies in the the
two-magnon continuum, potentially resulting in extra threshold singularities
in the magnon spectra.

###### Acknowledgements.

We would like to thank Pavel Maksimov for a prior collaboration on the earlier
attempt on this problem and for an important discussion concerning Kitaev-like
bond-dependent terms in 1D that has led us to expand on the general
anisotropic-exchange model for the zigzag chain. We are indebted to Leonie
Woodland and Radu Coldea for numerous conversations on the phenomenological
constrains for CoNb2O6, their implementation, and parameters of the model, as
well as for sharing their experimental results and for indispensable comments
and useful insights. We would like to thank Jeff Rau and Izabella Lovas for
helpful conversations. We are grateful to Shengtao Jiang for important
guidance regarding DMRG. This entire work, from conception to development,
execution, and writing, was supported by the U.S. Department of Energy, Office
of Science, Basic Energy Sciences under Award No. DE-SC0021221. We would like
to thank Aspen Center for Physics (A. L. C.) and KITP (A. L. C. and C. A. G.),
where parts of this work were completed. Aspen Center for Physics is supported
by National Science Foundation grant PHY-2210452. KITP is supported by the
National Science Foundation under Grants No. NSF PHY-1748958 and PHY-2309135.

##  Appendix A Model

###  A.1 Different parametrizations

For the two parametrizations of the exchange matrix in the crystallographic
reference frame $\\{a,b,c\\}$ in Fig. 2(a), the original one in Eq. (2) and
the one in the “ice-like” language in Eq. (5), the relations between exchanges
are given by

$\displaystyle J=\frac{J_{bb}+J_{cc}}{2}$ $\displaystyle,\
\Delta=\frac{2J_{aa}}{J_{bb}+J_{cc}},$ $\displaystyle
J_{z\pm}=\frac{J_{ac}}{\cos\varphi_{\alpha}}$ $\displaystyle,\
J_{\pm\pm}=\frac{J_{bb}-J_{cc}}{4\cos\varphi_{\alpha}},$
$\displaystyle\lambda_{z}=-\frac{(-1)^{\alpha}J_{ab}}{J_{ac}\tan\varphi_{\alpha}}$
$\displaystyle,\
\lambda_{\pm}=\frac{2(-1)^{\alpha}J_{bc}}{(J_{bb}-J_{cc})\tan\varphi_{\alpha}}.$
(51)

The matrix in the laboratory $\\{x_{0},y_{0},z_{0}\\}$ frame is obtained by
rotating the exchange matrix in the crystallographic frame by $\gamma$ about
$b$; see Fig. 2(b), using the rotation matrix

$\displaystyle\hat{\textbf{R}}_{\gamma}=\begin{pmatrix}\cos\gamma&0&-\sin\gamma\\\
0&1&0\\\ \sin\gamma&0&\cos\gamma\\\ \end{pmatrix}.$ (52)

The explicit relation of the matrix elements of the exchange matrix in the
laboratory frame to that in the crystallographic frame is

$\displaystyle J_{x_{0}x_{0}}$ $\displaystyle=J_{aa}\cos^{2}\\!\gamma-
J_{ac}\sin{2\gamma}+J_{cc}\sin^{2}\\!\gamma,$ $\displaystyle J_{y_{0}y_{0}}$
$\displaystyle=J_{bb},$ $\displaystyle J_{z_{0}z_{0}}$
$\displaystyle=J_{aa}\sin^{2}\\!\gamma+J_{ac}\sin{2\gamma}+J_{cc}\cos^{2}\\!\gamma,$
$\displaystyle J_{x_{0}y_{0}}$
$\displaystyle=J_{ab}\cos{\gamma}-J_{bc}\sin{\gamma},$ $\displaystyle
J_{x_{0}z_{0}}$
$\displaystyle=J_{ac}\cos{2\gamma}-\frac{1}{2}(J_{cc}-J_{aa})\sin{2\gamma},$
$\displaystyle J_{y_{0}z_{0}}$
$\displaystyle=J_{bc}\cos{\gamma}+J_{ab}\sin{\gamma}.$ (53)

The best-fit parameters in Ref. [44] can be translated to the exchanges in the
laboratory $\\{x_{0},y_{0},z_{0}\\}$ frame, see Eq. (12) and discussion in
Secs. II.4.1 and II.5.

Using Eqs. (53) and (12), we obtain exchanges in the crystallographic frame

$\displaystyle J_{aa}$
$\displaystyle=J_{x_{0}x_{0}}\cos^{2}\\!\gamma+J_{z_{0}z_{0}}\sin^{2}\\!\gamma=-1.05(1)\
\text{meV},$ $\displaystyle J_{bb}$ $\displaystyle=J_{y_{0}y_{0}}=-0.67(1)\
\text{meV},$ $\displaystyle J_{cc}$
$\displaystyle=J_{x_{0}x_{0}}\sin^{2}\\!\gamma+J_{z_{0}z_{0}}\cos^{2}\\!\gamma=-2.00(2)\
\text{meV},$ $\displaystyle J_{ab}$
$\displaystyle=J_{y_{0}z_{0}}\sin\gamma=-0.28(1)\ \text{meV},$ $\displaystyle
J_{ac}$ $\displaystyle=\frac{1}{2}(J_{z_{0}z_{0}}-J_{x_{0}x_{0}})\sin
2\gamma=-0.83(1)\ \text{meV},$ $\displaystyle J_{bc}$
$\displaystyle=J_{y_{0}z_{0}}\cos\gamma=-0.49(1)\ \text{meV}.$ (54)

Finally, combining Eqs. (51) and (54), gives parameters in the “ice-like”
parametrization in Eq. (13).

###  A.2 Out-of-plane angle

Real-space perturbation theory (RSPT) [64, 65, 66] allows to access the
effects of quantum fluctuations by expanding around the classical ground state
of the ferromagnetic Ising chain in various spin-flip processes. To avoid
unnecessary secondary details, we consider a simplified nearest-neighbor
exchange matrix

$\hat{\textbf{J}}_{\alpha}\\!=\\!\begin{pmatrix}0&(-1)^{\alpha}J_{x_{0}y_{0}}&0\\\
(-1)^{\alpha}J_{x_{0}y_{0}}&0&(-1)^{\alpha}J_{y_{0}z_{0}}\\\
0&(-1)^{\alpha}J_{y_{0}z_{0}}&{J}_{z_{0}z_{0}}\\\ \end{pmatrix},$ (55)

where $J_{z_{0}z_{0}}$ is the leading ferromagnetic Ising exchange and the two
staggered terms, $J_{x_{0}y_{0}}$ and $J_{y_{0}z_{0}}$, are perturbations. The
Hamiltonian can be written in terms of the spin ladder operators as follows

$\displaystyle\hat{\mathcal{H}}_{0}$
$\displaystyle=J_{z_{0}z_{0}}\sum_{\langle
ij\rangle}S_{i}^{z_{0}}S_{j}^{z_{0}},$ $\displaystyle\hat{V}_{xy}\\!$
$\displaystyle=\\!\frac{iJ_{x_{0}y_{0}}}{2}\sum_{\langle
ij\rangle}(-1)^{\alpha}\Big{\\{}S_{i}^{-}S_{j}^{-}-S_{i}^{+}S_{j}^{+}\Big{\\}},$
(56) $\displaystyle\hat{V}_{yz}\\!$
$\displaystyle=\\!\frac{iJ_{y_{0}z_{0}}}{2}\\!\sum_{\langle
ij\rangle}(-1)^{\alpha}\Big{\\{}\\!\left(S_{i}^{-}\\!-\\!S_{i}^{+}\\!\right)\\!S^{z_{0}}_{j}\\!+\\!\left(S_{j}^{-}\\!-\\!S_{j}^{+}\right)\\!S^{z_{0}}_{i}\\!\Big{\\}},$

where perturbations $\hat{V}_{xy}$ and $\hat{V}_{yz}$ generate double spin
flips and single spin flips, respectively. The ground state of the unperturbed
Hamiltonian in Eq. (56) is the ferromagnetic state $|0\rangle$ and its excited
states $|n\rangle$ are the states with $n$ spin flips.

Since we are interested in the deviations of the ordered moment from the Ising
axis, the lowest-order processes that induce a single-spin-flip state
$|1\rangle$ are in question. Notably, the single-spin-flip term acting on the
ground state vanishes identically because of its staggered form,
$\hat{V}_{yz}|0\rangle\\!=\\!0$, providing no spin tilt along the $y_{0}$
axis. The lowest non-zero contribution that yields the single-spin-flip state
is given by the second-order process involving both single- and double-spin-
flip terms in (56)

$|0\rangle\xrightarrow{\hat{V}_{xy}}|2\rangle\xrightarrow{\hat{V}_{yz}}|1\rangle,$
(57)

in which their mutually-canceling staggered form is important. Then, the
fluctuating ground state due to the process in (57) is

$\displaystyle|\widetilde{0}\rangle=|0\rangle-\frac{J_{x_{0}y_{0}}J_{y_{0}z_{0}}\sqrt{2S}}{8J_{z_{0}z_{0}}^{2}S\left(1-1/4S\right)}\,|1\rangle,$
(58)

which yields the angle of the spin tilt out of the $y_{0}z_{0}$ plane along
the $x_{0}$ axis, $\delta\gamma\\!\approx\\!\langle
S^{x_{0}}_{i}\rangle/\langle S^{z_{0}}_{i}\rangle$, for any site $i$

$\delta\gamma=-\frac{J_{x_{0}y_{0}}J_{y_{0}z_{0}}}{4J_{z_{0}z_{0}}^{2}S\left(1-1/4S\right)},$
(59)

where we used
$S^{x_{0}}_{i}\\!=\\!\left(S_{i}^{+}\\!+\\!S_{i}^{-}\right)\\!/2$ and
neglected higher-order corrections to the ground state from the two
consecutive double-spin-flips. The results in (58) and (59) are obtained for
the model (56) with arbitrary spin $S$, keeping higher-order $1/S$ terms such
as $(1-1/4S)$ factors, originating from the interaction of the nearest-
neighbor spin flips. The $S\\!=\\!1/2$ limit of (59) is listed in Eq. (9).

Importantly, the fluctuating ground state in (58) continues to respect the
glide symmetry and produces no tilt in the $y_{0}$ direction, $\langle
S^{y_{0}}_{i}\rangle\\!=\\!0$.

Figure 13: (a) The spin tilt angle by RSPT (59) and DMRG vs
$J_{x_{0}y_{0}}/|J_{z_{0}z_{0}}|$ for the simplified and full models for the
best-fit parameters for CoNb2O6. (b) The ratio of the DMRG and RSPT tilt
angles for the simplified model for two choices of $J_{y_{0}z_{0}}$, see the
text.

Given the analysis leading to the second-order RSPT result in (59), one can
argue that in order to produce a spin tilt, both $J_{x_{0}y_{0}}$ and
$J_{y_{0}z_{0}}$ terms are necessary in an arbitrary order of the theory,
because they have to cancel their staggered form. The higher-order corrections
to (59) also need to carry odd powers of each of the staggered terms because
they generate different number of the spin flips. These considerations can be
expected to remain valid for a general model in Eq. (6), which contains other
non-staggered spin-flip terms.

As is discussed in Sec. II.4.1 and Sec. II.4.2 and is clear from Eq. (59), the
spin tilt angle vanishes for the choice of $J_{x_{0}y_{0}}\\!=\\!0$ made in
Refs. [43, 44], rendering quantum corrections to the classical Ising spin
direction zero. One can verify the accuracy of the second-order perturbative
result for the tilt angle in (59) and elucidate the role of the higher-order
fluctuations with the help of the unbiased DMRG calculations for the ground
state.

The DMRG simulations were performed in chains of up to $5000$ sites using the
ITensor library [67]. We have employed two strategies, the “scan” with a
slowly varying $J_{x_{0}y_{0}}$ along the chain, which provides a real-space
variation of the tilt angle, and “no-scan,” in which the tilt angle is
measured in the middle of the chain away from the edges for each individual
$J_{x_{0}y_{0}}$ value. Both approaches yield numerically indistinguishable
results. Because of the Ising nature of the model, a very good convergence is
reached with low bond dimensions and small number of the DMRG sweeps [84].

In Figure 13(a), we show the RSPT result (59) for $S\\!=\\!1/2$ and for
$J_{z_{0}z_{0}}$ and $J_{y_{0}z_{0}}$ from Eq. (12), which correspond to the
best-fit parameters for CoNb2O6, as a function of
$J_{x_{0}y_{0}}/|J_{z_{0}z_{0}}|$. It is shown together with the DMRG results
for the same $S\\!=\\!1/2$ simplified model (55), for which Eq. (59) was
derived. In addition to that, DMRG results for the full model of CoNb2O6 in
Eqs. (7) and (14) for the best-fit parameters from Eqs. (12) and (15) in Secs.
II.5 are shown.

Clearly, the numerical results for $\delta\gamma$ vanish identically for
$J_{x_{0}y_{0}}\\!=\\!0$, in accord with the discussion above. For the
simplified model, the agreement of the slope of the DMRG tilt angle with that
of the RSPT is close, but not precise, as is also demonstrated in Fig. 13(b),
which shows the ratio of the angles. The majority of this difference can be
attributed to the next-order correction, corresponding to the forth-order
process

$|0\rangle\xrightarrow{\hat{V}_{xy}}|2\rangle\xrightarrow{\hat{V}_{yz}}|3\rangle\xrightarrow{\hat{V}_{yz}}|2\rangle\xrightarrow{\hat{V}_{yz}}|1\rangle,$
(60)

which is $\mathcal{O}(J_{x_{0}y_{0}}J_{y_{0}z_{0}}^{3})$, conforming to the
rules proposed for the higher-order corrections that are discussed above.

One can verify the $J_{y_{0}z_{0}}^{3}$-order of the correction to the slope
by changing the numerical value of $J_{y_{0}z_{0}}$ as is shown in Fig. 13(b).
Here, the reduction of $J_{y_{0}z_{0}}$ by $1/2$ leads to an eightfold
decrease of the correction in the DMRG result, all in a close accord with the
expectations of the odd powers of each staggered term outlined above.

The DMRG results for the best-fit parameters in the full model of CoNb2O6 in
Fig. 13(a) show a qualitative agreement with the perturbative result (59) for
the simplified model (55), but also a substantial quantitative difference. It
can be attributed to the fluctuations induced by the other spin-flip terms
present in the full model.

Lastly, as is discussed in Sec. II.4.2, the phenomenological constraint on the
spin direction in the model (6) allows the “residual” $J_{x_{0}y_{0}}$ and
$J_{x_{0}z_{0}}$ terms to be present, but exactly compensating each others’
spin tilting, leaving the physical Ising direction intact. The perturbative
result in Eq. (59) together with the classical energy minimization result in
Eq. (8) immediately suggest an explicit connection between the two terms:
$J_{x_{0}z_{0}}\\!=\\!-\delta\gamma J_{z_{0}z_{0}}$, see Sec. II.4.2 for more
detail. The DMRG results shown in Fig. 4 in that Section, demonstrating a
relation between $J_{x_{0}y_{0}}$ and $J_{x_{0}z_{0}}$, were calculated using
the best-fit parameters of the full model of CoNb2O6. For each fixed value of
$J_{x_{0}y_{0}}$, we performed a DMRG scan vs $J_{x_{0}z_{0}}$ to identify the
value of $J_{x_{0}z_{0}}$ that corresponds to an exact compensation of the
spin tilt away from the Ising axis.

##  Appendix B Spin-wave theory

We consider a standard Holstein-Primakoff (HP) spin representation with the
quantization axis $z$

$S_{i}^{z}=S-n_{i},\ \ \ \ S_{i}^{+}=a_{i}\sqrt{2S-n_{i}},$ (61)

where $a_{i}(a_{i}^{\dagger})$ are bosonic operators,
$n_{i}\\!=\\!a_{i}^{\dagger}a_{i}$, and $i$ is the site index. The expansion
of the square roots in (61) in powers of $n_{i}/2S$ results in the bosonic
Hamiltonian

$\mathcal{\hat{H}}=\mathcal{\hat{H}}^{(0)}+\mathcal{\hat{H}}^{(1)}+\mathcal{\hat{H}}^{(2)}+\mathcal{\hat{H}}^{(3)}+\mathcal{\hat{H}}^{(4)}+\mathcal{O}(S^{-1}),$
(62)

where the $n$th term $\mathcal{\hat{H}}^{(n)}$ contains the $n$th power of the
bosonic operators and carries an explicit $S^{2-n/2}$ factor, constituting a
$1/S$-expansion for a given problem.

The first term $\hat{\mathcal{H}}^{(0)}$ in such an expansion is the classical
energy, and $\hat{\mathcal{H}}^{(1)}$ should vanish upon the classical energy
minimization. The quadratic term $\hat{\mathcal{H}}^{(2)}$ is the harmonic
part of the expansion that yields the LSWT and magnon energy spectrum. The
lowest $1/S$-order corrections to the LSWT originate from the two anharmonic
terms in the expansion, $\hat{\mathcal{H}}^{(3)}$ and
$\hat{\mathcal{H}}^{(4)}$, describing three- and four-magnon interaction,
respectively.

###  B.1 Classical energy

The classical energy of the field-polarized paramagnetic phase considered in
Sec. III.1 is easily obtained from (19) and is given by
$E_{\text{cl}}/N\\!=\\!-SH+S^{2}(J_{y_{0}y_{0}}+J_{2})$, with the linear
$\hat{\mathcal{H}}^{(1)}$ from (20) vanishing because of the staggered
structure of the bond-dependent terms. This is a common situation for
collinear states that do not require energy minimization and cannot indicate
their phase boundaries from the classical consideration alone.

One standard approach is to proceed directly with the harmonic term
$\hat{\mathcal{H}}^{(2)}$, develop LSWT as in Sec. III.2, and obtain the value
of the critical field from the condition of stability of the magnon spectrum
in (26).

The other approach is to consider the ordered phase of CoNb2O6 in the
transverse field in the classical limit and find the critical field of a
transition to the fully polarized state from the minimization of its energy.
In such a state spins are tilted away from the field toward the Ising axis in
the $y_{0}z_{0}$ ($zx$) plane; see Fig. 2(c). Denoting this angle as $\theta$,
straightforward algebra in (19) yields

$\displaystyle\frac{E_{\text{cl}}}{N}=$ $\displaystyle-
SH\cos\theta+S^{2}(J_{y_{0}y_{0}}+J_{2})$
$\displaystyle-S^{2}(J_{y_{0}y_{0}}-J_{z_{0}z_{0}}+J_{2}-J_{2z_{0}})\sin^{2}\theta.$
(63)

Minimizing it with respect to $\theta$ gives $\cos\theta\\!=\\!H/H_{c}$ with
the critical field given in Eq. (17).

###  B.2 Linear spin-wave theory

The LSWT is based on the lowest-order expansion in (61) in the Hamiltonian
(19), leading to

$\displaystyle\hat{\mathcal{H}}^{(2)}$
$\displaystyle=\frac{S}{2}\sum_{i}\Big{\\{}2\big{(}H/S-2(J_{y_{0}y_{0}}+J_{2})\big{)}n_{i}$
$\displaystyle+\Big{(}(J_{z_{0}z_{0}}+J_{x_{0}x_{0}})a_{i}^{\dagger}a_{i+1}^{\phantom{{\dagger}}}+(J_{2z}+J_{2})a_{i}^{\dagger}a_{i+2}^{\phantom{{\dagger}}}$
(64)
$\displaystyle+(J_{z_{0}z_{0}}\\!-\\!J_{x_{0}x_{0}})a_{i}^{\dagger}a_{i+1}^{\dagger}+(J_{2z}-J_{2})a_{i}^{\dagger}a_{i+2}^{\dagger}+\text{H.c.}\Big{)}\Big{\\}},$

Using Fourier transformation (21) gives the harmonic Hamiltonian in the
canonical form (22). The standard Bogolyubov transformation diagonalizes it
with the $u_{\bf k}$ and $v_{\bf k}$ parameters given explicitly as

$\displaystyle u_{\bf k}=\sqrt{\frac{A_{\bf k}+\varepsilon_{\bf
k}}{2\varepsilon_{\bf k}}},\ v_{\bf k}=\text{sgn}(B_{\bf k})\sqrt{\frac{A_{\bf
k}-\varepsilon_{\bf k}}{2\varepsilon_{\bf k}}},$ (65)

with $A_{\bf k}$ and $B_{\bf k}$ from Eq. (23) and $\varepsilon_{\bf k}$ from
Eq. (25). The resultant diagonal form of the LSWT Hamiltonian is

$\mathcal{\hat{H}}^{(2)}=\sum_{\bf k}\Big{\\{}\varepsilon_{\bf k}b_{\bf
k}^{\dagger}b_{\bf k}^{\phantom{{\dagger}}}+\frac{1}{2}\left(\varepsilon_{\bf
k}-A_{\bf k}\right)\Big{\\}},$ (66)

where the magnon energy in (25) can also be written as

$\displaystyle\varepsilon_{\bf k}^{2}\\!=$
$\displaystyle\Big{[}H\\!-\\!2S\left(J_{y_{0}y_{0}}\\!-\\!J_{x_{0}x_{0}}\gamma^{(1)}_{\bf
k}+J_{2}\left(1-\gamma^{(2)}_{\bf k}\right)\right)\Big{]}$
$\displaystyle\times$
$\displaystyle\Big{[}H\\!-\\!2S\left(J_{y_{0}y_{0}}\\!+\\!J_{2}\\!-\\!J_{z_{0}z_{0}}\gamma_{\bf
k}^{(1)}\\!-\\!J_{2z_{0}}\gamma_{\bf k}^{(2)}\right)\Big{]}.$ (67)

###  B.3 Non-linear spin-wave theory

####  B.3.1 Cubic terms

Using the leading-order HP expansion in the odd part of the Hamiltonian (20)
leads to the cubic term

$\mathcal{\hat{H}}^{(3)}\\!=\\!-J_{y_{0}z_{0}}\sqrt{\frac{S}{2}}\sum_{i}\\!\Big{(}(-1)^{i}n_{i}\big{(}a^{\dagger}_{i+1}\\!-a^{\dagger}_{i-1}\big{)}\\!+\text{H.c.}\\!\Big{)}.$
(68)

The Fourier transformation (21) in Eq. (68) yields

$\displaystyle\mathcal{\hat{H}}^{(3)}=J_{y_{0}z_{0}}\sqrt{\frac{2S}{N}}\sum_{{\bf
k},{\bf q}}\left(\bar{\gamma}_{\bf p}a^{\dagger}_{\bf p}a^{\dagger}_{\bf
q}a_{\bf k}^{\phantom{{\dagger}}}+\text{H.c.}\right),$ (69)

with ${\bf p}={\bf k}-{\bf q}+{\bf Q}$ and $\bar{\gamma}_{\bf p}=i\sin(p{\sf
c}_{0})$. The Bogolyubov transformation with symmetrization give

$\displaystyle\mathcal{\hat{H}}^{(3)}$
$\displaystyle=\frac{1}{2!\sqrt{N}}\sum_{{\bf q}+{\bf k}+{\bf p}={\bf
Q}}\Big{(}\Phi_{{\bf qk};{\bf p}}b^{\dagger}_{\bf q}b^{\dagger}_{\bf
k}b_{-{\bf p}}^{\phantom{{\dagger}}}+\text{H.c.}\Big{)}$
$\displaystyle+\frac{1}{3!\sqrt{N}}\sum_{{\bf q}+{\bf k}+{\bf p}={\bf
Q}}\Big{(}\Xi_{{\bf qkp}}b^{\dagger}_{\bf q}b^{\dagger}_{\bf
k}b^{\dagger}_{\bf p}+\text{H.c.}\Big{)},$ (70)

with the decay and source vertices, $\Phi_{{\bf qk};{\bf p}}$ and $\Xi_{{\bf
qkp}}$,

$\displaystyle\Phi_{{\bf qk};{\bf
p}}\\!=\\!\sqrt{2S}\,J_{y_{0}z_{0}}\widetilde{\Phi}_{{\bf qk};{\bf p}},\ \
\Xi_{{\bf qk}{\bf p}}\\!=\\!\sqrt{2S}\,J_{y_{0}z_{0}}\widetilde{\Xi}_{{\bf
qk}{\bf p}},$ (71)

and the dimensionless vertices given by

$\displaystyle\widetilde{\Phi}_{{\bf qk};{\bf p}}$
$\displaystyle=\bar{\gamma}_{\bf k}(u_{\bf k}+v_{\bf k})(u_{\bf q}u_{\bf
p}+v_{\bf q}v_{\bf p})$ $\displaystyle+\bar{\gamma}_{\bf q}(u_{\bf q}+v_{\bf
q})(u_{\bf k}u_{\bf p}+v_{\bf k}v_{\bf p})$ $\displaystyle+\bar{\gamma}_{\bf
p}(u_{\bf p}+v_{\bf p})(u_{\bf k}v_{\bf q}+v_{\bf k}u_{\bf q}),$ (72)
$\displaystyle\widetilde{\Xi}_{{\bf qk}{\bf p}}$
$\displaystyle=\bar{\gamma}_{\bf k}(u_{\bf k}+v_{\bf k})(u_{\bf q}v_{\bf
p}+v_{\bf q}u_{\bf p})$ $\displaystyle+\bar{\gamma}_{\bf q}(u_{\bf q}+v_{\bf
q})(u_{\bf k}v_{\bf p}+v_{\bf k}u_{\bf p})$ $\displaystyle+\bar{\gamma}_{\bf
p}(u_{\bf p}+v_{\bf p})(u_{\bf k}v_{\bf q}+v_{\bf k}u_{\bf q}).$ (73)

The decay and source vertices in (70) are umklapp-like, with the momentum
conserved up to the ${\bf Q}$-vector.

The resulting lowest-order self-energies are

$\displaystyle\Sigma^{(d)}({\bf k},\omega)\\!$
$\displaystyle=\\!\frac{1}{2N}\sum_{\bf q}\frac{|\Phi_{{\bf q},{\bf k}-{\bf
q}+{\bf Q};-{\bf k}}|^{2}}{\omega-\varepsilon_{\bf q}-\varepsilon_{{\bf
k}-{\bf q}+{\bf Q}}+i0^{+}},$ (74) $\displaystyle\Sigma^{(s)}({\bf
k},\omega)\\!$ $\displaystyle=\\!-\frac{1}{2N}\sum_{\bf q}\frac{|\Xi_{{\bf
q},-{\bf k}-{\bf q}+{\bf Q},{\bf k}}|^{2}}{\omega+\varepsilon_{\bf
q}+\varepsilon_{-{\bf k}-{\bf q}+{\bf Q}}-i0^{+}}.$ (75)

####  B.3.2 Quartic terms

The four-boson terms are obtained from the higher-order expansion in the HP
transformation (61)

$S_{i}^{+}\approx\sqrt{2S}\left(a_{i}-\frac{n_{i}a_{i}}{4S}\right),\hskip
14.22636ptS^{z}_{i}=S-n_{i}.$ (76)

Using this expansion (76), the quartic terms come from

$\displaystyle S_{i}^{x}S_{j}^{x}$
$\displaystyle\rightarrow-\frac{1}{8}\Big{(}(a_{i}^{\dagger}+a_{i})n_{j}a_{j}+(i\leftrightarrow
j)+\text{H.c.}\Big{)},$ $\displaystyle S_{i}^{y}S_{j}^{y}$
$\displaystyle\rightarrow-\frac{1}{8}\Big{(}(a_{i}^{\dagger}-a_{i})n_{j}a_{j}+(i\leftrightarrow
j)+\text{H.c.}\Big{)},$ $\displaystyle S_{i}^{z}S_{j}^{z}$
$\displaystyle\rightarrow n_{i}n_{j}.$ (77)

The decoupling of them uses the real-space HF averages

$\displaystyle n$ $\displaystyle=\langle
a_{i}^{\dagger}a_{i}^{\phantom{{\dagger}}}\rangle\\!=\\!\sum_{\bf k}v_{\bf
k}^{2},\hskip 5.69046ptm_{n}\\!=\\!\langle
a_{i}^{\dagger}a_{j}^{\phantom{{\dagger}}}\rangle\\!=\\!\sum_{\bf
k}\gamma_{\bf k}^{(n)}v_{\bf k}^{2},$ $\displaystyle\delta$
$\displaystyle=\langle a_{i}^{2}\rangle\\!=\\!\sum_{\bf k}u_{\bf k}v_{\bf
k},\hskip 5.69046pt\Delta_{n}\\!=\\!\langle a_{i}a_{j}\rangle\\!=\\!\sum_{\bf
k}\gamma_{\bf k}^{(n)}u_{\bf k}v_{\bf k},$ (78)

with the index $n\\!=\\!1(2)$ for the nearest and next-nearest neighbors,
yielding

$\displaystyle S^{x}_{i}S^{x}_{j}\left(S^{y}_{i}S^{y}_{j}\right)\rightarrow$
$\displaystyle-\frac{1}{2}\biggl{[}(m_{n}\pm\Delta_{n})(n_{i}+n_{j})$
$\displaystyle+\frac{1}{4}(\Delta_{n}\pm
m_{n})\left(a_{i}^{\dagger}a_{i}^{\dagger}+a_{j}^{\dagger}a_{j}^{\dagger}+\text{H.c.}\right)$
$\displaystyle+\left(n\pm\frac{\delta}{2}\right)\left(a_{i}^{\dagger}a_{j}^{\phantom{{\dagger}}}+\text{H.c.}\right)$
$\displaystyle+\left(\frac{\delta}{2}\pm
n\right)\left(a_{i}^{\dagger}a_{j}^{\dagger}+\text{H.c.}\right)\biggl{]},$
(79) $\displaystyle S_{i}^{z}S_{j}^{z}\\!\rightarrow
n(n_{i}+n_{j})+\\!\big{(}m_{n}a_{i}^{\dagger}a_{j}+\Delta_{n}a_{i}^{\dagger}a_{j}^{\dagger}+\text{H.c.}\big{)}.$
(80)

Using (79) and (80) in the quartic Hamiltonian from the even part of the model
(19), followed by the Fourier transformation, gives a correction to the LSWT
model

$\displaystyle\delta\hat{\mathcal{H}}^{(4)}=\sum_{\bf k}\Big{\\{}\delta A_{\bf
k}a_{\bf k}^{\dagger}a_{\bf k}^{\phantom{{\dagger}}}-\frac{\delta B_{\bf
k}}{2}\left(a_{\bf k}^{\dagger}a_{-{\bf
k}}^{\dagger}+\text{H.c}\right)\Big{\\}},\quad\quad$ (81)

with

$\displaystyle\delta A_{\bf k}$
$\displaystyle=-J_{x_{0}x_{0}}\Big{(}m_{1}-\Delta_{1}+\left(n-\delta/2\right)\gamma_{\bf
k}^{(1)}\Big{)}$ $\displaystyle+2J_{y_{0}y_{0}}\Big{(}n+m_{1}\gamma_{\bf
k}^{(1)}\Big{)}$ $\displaystyle-
J_{z_{0}z_{0}}\Big{(}m_{1}+\Delta_{1}+\left(n+\delta/2\right)\gamma_{\bf
k}^{(1)}\Big{)}$ $\displaystyle-
J_{2}\Big{(}m_{2}-\Delta_{2}-2n+\left(n-\delta/2-2m_{2}\right)\gamma_{\bf
k}^{(2)}\Big{)}$ $\displaystyle-
J_{2z}\Big{(}m_{2}+\Delta_{2}+\left(n+\delta/2\right)\gamma_{\bf
k}^{(2)}\Big{)},$ (82) $\displaystyle\delta B_{\bf k}$
$\displaystyle=-J_{x_{0}x_{0}}\biggl{(}\frac{1}{2}(m_{1}-\Delta_{1})+(n-\delta/2)\gamma_{\bf
k}^{(1)}\biggl{)}$ $\displaystyle-2J_{y_{0}y_{0}}\Delta_{1}\gamma_{\bf
k}^{(1)}$
$\displaystyle+J_{z_{0}z_{0}}\biggl{(}\frac{1}{2}(m_{1}+\Delta_{1})+\left(n+\delta/2\right)\gamma_{\bf
k}^{(1)}\biggl{)}$ $\displaystyle-
J_{2}\biggl{(}\frac{1}{2}(m_{2}-\Delta_{2})+(n-\delta/2+2\Delta_{2})\gamma_{\bf
k}^{(2)}\biggl{)}$
$\displaystyle+J_{2z_{0}}\biggl{(}\frac{1}{2}(m_{2}+\Delta_{2})+\left(n+\delta/2\right)\gamma_{\bf
k}^{(2)}\biggl{)}.$ (83)

Finally, using Bogolyubov transformation in (81), yields the
$\omega$-independent $1/S$ energy correction in (29).

###  B.4 Self-consistent Hartree-Fock method

Here we provide some further details on the self-consistent HF method
discussed in Sec. III.4.

The set of the real-space HF averages in (78),
$\\{\mbox{HFs}\\}\\!=\\!\\{n,\delta,m_{1},m_{2},\Delta_{1},\Delta_{2}\\}$, are
found iteratively for each fixed field value, starting from
$H_{0}\\!>\\!H_{c}$ and proceeding by decreasing the field with a small step
$\Delta H$.

The initial set of the HF averages, $\\{\mbox{HFs}\\}_{n}^{0}$, for a field
$H_{n}$ in such a sequence of fields is taken from the final (converged) set
of the HF averages from the previous field value, $\\{\mbox{HFs}\\}_{n-1}^{\rm
final}$. The exception is the very first field $H_{0}$, which we choose large
enough for the LSWT averages from (78) to be a good starting point, so we use
$\\{\mbox{HFs}\\}_{0}^{0}\\!=\\!\\{\mbox{HFs}\\}_{0}^{\rm LSWT}$ for it.

For any field $H_{n}$, the self-consistent iterations follow the cycle shown
in (33), with the steps for each subsequent iteration summarized as follows:

  * 1)

At each $i$th step, the $\\{\mbox{HFs}\\}_{n}^{i-1}$ averages give the
quartic-term contributions to the harmonic theory, $\delta\bar{A}_{\bf k}$ and
$\delta\bar{B}_{\bf k}$, according to the expressions in (82) and (83). The
LSWT-like eigenvalue problem of the same form as in Eq. (22) with
$\bar{A}_{\bf k}\\!=\\!A_{\bf k}+\delta\bar{A}_{\bf k}$ and $\bar{B}_{\bf
k}\\!=\\!B_{\bf k}+\delta\bar{B}_{\bf k}$ yields the new set of the Bogolyubov
parameters $\bar{u}_{\bf k}$ and $\bar{v}_{\bf k}$.

  * 2)

Using $\bar{u}_{\bf k}$ and $\bar{v}_{\bf k}$, a new (temporary) set of HF
averages, $\\{\mbox{HFs}\\}_{n}^{i_{t}}$, is calculated using Eq. (78).

  * 3)

The input of the HF averages for the next iteration is updated using
$\\{\mbox{HFs}\\}_{n}^{i}=\alpha\\{\mbox{HFs}\\}_{n}^{i_{t}}+(1-\alpha)\\{\mbox{HFs}\\}_{n}^{i-1}$,
where $\alpha\\!\ll\\!1$ ensures a smooth convergence.

  * 4)

The cycle of the steps from 1) to 3) is continued until a numerical
convergence in the HF averages is reached by meeting a tolerance $\epsilon$
between the two subsequent iterations. At this step, the final set
$\\{\mbox{HFs}\\}_{n}^{\rm final}$ for the field $H_{n}$ is defined.
Obviously, it also yields the SCHF magnon eigenenergies
$\bar{\varepsilon}_{\bf k}$ and Bogolyubov parameters $\bar{u}_{\bf k}$ and
$\bar{v}_{\bf k}$ used in our results in Sec. III.4 and Sec. IV.

  * 5)

For the next field $H_{n+1}$, the cycle starts at the step 1) with the
converged set of $\\{\mbox{HFs}\\}_{n}^{\rm final}$ used as the initial
condition.

In this work, we have used $H_{0}\\!=\\!3H_{c}$, $\Delta
H\\!=\\!3\\!\times\\!10^{-3}\ \text{meV}$, $\alpha\\!=\\!0.01$, and
$\epsilon\\!=\\!10^{-7}$. The stability of this procedure was verified by
varying all of these parameters to ensure independence of the results.

Figure 14: Representative Hartree-Fock averages vs $H$, LSWT (dashed lines)
and SCHF (solid lines), respectively.

Our Fig. 14 illustrates that the logarithmically divergent behavior of the HF
averages gets regularized using the SCHF method and that their calculation is
successfully extended below the LSWT value of $H_{c}$.

##  Appendix C Results

###  C.1 DMRG calculations of the gap

The DMRG simulations for the lowest energy gap in the excitation spectrum of
the best-fit model of CoNb2O6 as a function of the field that is shown in Fig.
8 in Sec. IV were performed in the chains of up to $500$ sites. The gap was
obtained by calculating the ground state and the first excited state using the
ITensor library [67] for each value of the field. While a good convergence was
easily reached for the ground state [84], a larger number of sweeps (about
$75$) were needed to reach the same accuracy when computing the first excited
state.

Figure 15: (a) Energy gap vs field by DMRG for the chains of different length.
(b) The expectation value of the spin component $\langle S_{i}^{z_{0}}\rangle$
in the ground state along the long chain in a DMRG scan vs field. Vertical
line indicates the critical field. Inset shows the second derivative of
$\langle S_{i}^{z_{0}}\rangle$ vs field.

Fig. 15(a) shows the excitation gap in a narrow field region near the critical
field for the chains of lengths up to $L\\!=\\!1000$. One can see that the
finite-size effects are appreciable only near the minimum of the gap. Using
the $1/L$ extrapolation of the data in Fig. 15(a) for small gaps, we have
verified that the $L\\!\rightarrow\infty$ gap vanishes at about $4.52(1)$ T,
the value which is also consistent with the linear $|B-B_{c}|$ extrapolation
of the $L\\!=\\!1000$ data above the critical field. It is worth noting that
in the ordered phase, the finite-size effects lift the ground state
degeneracy, with the resultant splitting that can be confused with the actual
excitation gap close to the critical field, as both reach
$\mathcal{O}(10^{-3})$ meV for $L\\!=\\!1000$.

One can also corroborate these results for the critical field using the
ground-state DMRG calculations. Our Fig. 15(b) shows the DMRG scan in the
ground state with a slowly varying field in a chain of $5000$ sites. The
expectation value of the spin component in the Ising direction $\langle
S_{i}^{z_{0}}\rangle$ suggests a second order phase transition with a critical
field at about $4.527(1)\ \mathrm{T}$. This value is determined from the
inflection point of the curve, with the change on the sign of the second
derivative shown in the inset. This result is consistent with the value
suggested by the analysis of the gap from Fig. 15(a).

###  C.2 Two-magnon kinematics

Here we briefly discuss two aspects of the two-magnon kinematics in the
context of the 1D spin model of CoNb2O6: the structure of the two-magnon
continuum and kinematics of the magnon decay.

####  C.2.1 The two-magnon continuum

At any given momentum ${\bf k}$, the ${\bf Q}$-shifted two-magnon continuum is
defined within the energy range

$\displaystyle E_{{\bf k},{\bf q}^{*}}^{\text{min}}\leq E_{{\bf k},{\bf
q}}=\bar{\varepsilon}_{\bf q}+\bar{\varepsilon}_{{\bf k}-{\bf q}+{\bf Q}}\leq
E_{{\bf k},{\bf q}^{*}}^{\text{max}},$ (84)

where the boundaries $E_{{\bf k},{\bf q}^{*}}^{\text{min}}$ and $E_{{\bf
k},{\bf q}^{*}}^{\text{max}}$ should be found from the extremum condition
$\left(\partial E_{{\bf k},{\bf q}}/\partial{\bf q}\right)|_{{\bf
q}^{*}}\\!=\\!0$, which translates to the requirement on the group velocities
of the two magnons to be equal for the momentum ${\bf q}^{*}$

$\displaystyle\bar{\bf v}_{{\bf q}^{*}}=\bar{\bf v}_{{\bf k}-{\bf q}^{*}+{\bf
Q}}.$ (85)

Although, generally, such conditions may require a numerical solution, one
class of them, which often describes a majority of the extrema in the two-
magnon continua [46, 71], is easy to find analytically. The equivalence of the
magnon velocities in (85) is automatically satisfied when their momenta are
identical up to a set of the reciprocal lattice vectors. In our case, using
that the smallest reciprocal lattice vector is ${\bf G}\\!=\\!2{\bf Q}$, one
obtains two solutions, ${\bf q}^{*}_{\pm}\\!=\\!({\bf k}\pm{\bf Q})/2$,
referred to as the “equivalent magnon” solutions below. The energies of the
two-magnon continua for them are $E_{{\bf k},{\bf
q}^{*}_{\pm}}\\!=\\!2\bar{\varepsilon}_{{\bf q}^{*}_{\pm}}$.

Figure 16: (a) The two-magnon DoS for the best-fit model of CoNb2O6 and
$B\\!=\\!7$ T ($H\\!=\\!0.8H_{c}$) as in Fig. 9, focusing on the vicinity of
the ${\bf\Gamma}$ point; dashed lines are the equivalent-magnon solutions
$E_{{\bf k},{\bf q}^{*}_{\pm}}$. (b) The energy of the two-magnon continuum
$E_{{\bf k},{\bf q}}$ vs ${\bf q}$ for several ${\bf k}$, with the extrema
corresponding to different solutions for ${\bf q}^{*}$ indicated, see text.

Given the relative simplicity of the magnon dispersion in either the LSWT or
SCHF treatment of the model for CoNb2O6, and because of its 1D character,
which permits only maxima and minima, these two solutions indeed describe the
two edges of the two-magnon continuum and completely exhaust the number
singularities in it for most of the ${\bf k}$ values, as one can see in Fig.
9, Sec. IV. Our Fig. 16(a) shows the intensity plot of the two-magnon DoS for
the best-fit model of CoNb2O6 for $B\\!=\\!7$ T from Fig. 9, focusing instead
on the vicinity of the ${\bf\Gamma}$ point, which provides a clearer view of
the significantly richer structure of the continuum in that region.

While additional singularities in the continuum in the 2D and 3D cases are
common and are usually associated with the saddle points within it [46, 71],
the appearance of multiple singularities for the 1D model of CoNb2O6 is
somewhat of a surprise. In Fig. 16(a), one can see that the equivalent-magnon
extrema cease to be the absolute minima and maxima of the continuum and are
overtaken by two different ones in the proximity of the $\Gamma$ point.

Another insight into the structure of the continuum is offered by Fig. 16(b),
which shows the ${\bf q}$-cut of the two-magnon continuum for ${\bf
k}\\!=\\!0$ with the additional minima and maxima, demonstrating that the
equivalent-magnon extremum at ${\bf q}^{*}_{\pm}\\!=\\!\pm{\bf Q}/2$ is now a
local minimum.

Counterintuitively, this unusual complexity is an outcome of the simplicity
and high symmetry of the magnon spectrum. As is discussed in Sec. IV.2, away
from the small-gap regime near the critical field, the magnon energy can be
well-described by the nearest-neighbor hopping approximation,
$\bar{\varepsilon}_{\bf q}\\!\approx\\!{\sf E_{0}}+{\sf J_{1}}\gamma_{\bf
q}^{(1)}$, which naturally explains the bow-tie form of the continuum at the
${\bm{\Gamma}}$ point, as the nearest-neighbor hopping terms of
$\bar{\varepsilon}_{\bf q}$ and $\bar{\varepsilon}_{{\bf k}-{\bf q}+{\bf Q}}$
in the continuum (84) cancel precisely at ${\bf k}\\!=\\!0$.

However, the next-nearest-neighbor exchanges (24) and relativistic form of the
magnon dispersion (25) produce small, but essential further-neighbor hoppings.
Specifically, the finite width of the continuum at the ${\bm{\Gamma}}$ point
can only be provided by the effective even-neighbor hoppings. One can verify
that the features shown in Fig. 16 can be closely reproduced by the following
approximation for the magnon energy

$\displaystyle\bar{\varepsilon}_{\bf q}\approx{\sf E_{0}}+{\sf
J_{1}}\gamma_{\bf q}^{(1)}+{\sf J_{2}}\gamma_{\bf q}^{(2)}+{\sf
J_{4}}\gamma_{\bf q}^{(4)},$ (86)

using ${\sf E_{0}}\\!=\\!2.107$, ${\sf J_{1}}\\!=\\!-1.4$, ${\sf
J_{2}}\\!=\\!-0.0065$, and ${\sf J_{4}}\\!=\\!-0.01$, all in meV, where
$\gamma_{\bf k}^{(n)}\\!=\\!\cos(nk_{c}{\sf c}_{0})$ as before.

Using the form in (86), a simple algebra yields the maxima of the ${\bf
k}\\!=\\!0$ continuum in Fig. 16(b) at the momenta ${\bf
q}^{*,n}_{4\pm}\\!=\\!\pm\frac{1}{2}\arccos(-{\sf J_{2}}/4{\sf J_{4}})+n\pi$,
explicating the essential role of the further-neighbor terms in the additional
extrema of the continuum. For small ${\bf k}$, the left pair of these maxima
shifts in ${\bf q}$ and up in energy linearly with ${\bf k}$ and remains the
absolute maxima for a range of ${\bf k}$, while the right pair shifts down in
energy, becomes the local maxima, and then ceases to be extremal at the larger
${\bf k}$. In Fig. 16(a), they correspond to the upper singularity, which
merges with the one from the equivalent magnons, and to the one that enters
the continuum and annihilates with the nearly flat singularity, respectively.

The last and the most curious is the “flat” singularity, which is the absolute
minimum of the continuum in Fig. 16(a) at ${\bf k}\\!=\\!0$. The corresponding
minima of $E_{{\bf k},{\bf q}}$ in Fig. 16(b) can be found at ${\bf
q}^{*}_{0}\\!\approx\\!-b_{0}{\bf k}$ and ${\bf q}^{*}_{\bf
Q}\\!\approx\\!{\bf Q}+b_{0}{\bf k}$, where $b_{0}\\!=\\!\frac{1}{2}({\sf
J_{1}}/\widetilde{{\sf J}}_{{\sf 0}}\\!-\\!1)$ with $\widetilde{{\sf J}}_{{\sf
0}}\\!=\\!4({\sf J_{2}}+4{\sf J_{4}})$. At ${\bf k}\\!=\\!0$, one magnon in
$E_{0,{\bf q}^{*}}$ is at ${\bf q}^{*}_{0}\\!=\\!0$ and the other is at ${\bf
Q}$, precisely at the minimum and the maximum of the single-magnon band. This
arrangement guarantees that the velocities on both sides of Eq. (85) are zero,
fulfilling the extremum condition. It also explains the flatness of the
singularity in Fig. 16(a), as the continuum energy is $E_{{\bf k},{\bf
q}^{*}_{0}}\\!\approx\\!\bar{{\sf E}}_{\sf 0}+2\widetilde{{\sf J}}_{{\sf
0}}b_{0}^{2}{\bf k}^{2}$, with $\bar{{\sf E}}_{\sf 0}\\!=\\!2({\sf E_{0}}+{\sf
J_{2}}+{\sf J_{4}})$.

This completes the consideration of the richer set of the Van Hove
singularities in the two-magnon continuum near the ${\bf\Gamma}$ point in our
model.

####  C.2.2 The kinematics of the two-magnon decay

Figure 17: (a) Schematics of the two-magnon decay of a magnon for ${\bf
k}\\!>\\!{\bf k}^{*}$; solid line is $\bar{\varepsilon}_{\bf k}$ and shaded
area is the continuum. (b) The set of decay ${\bf q}$-points vs ${\bf k}$;
vertical lines are marking decay thresholds ${\bf k}^{*}$ and ${\bf G}-{\bf
k}^{*}$. Arrows indicate the shift by the reciprocal lattice vector $-{\bf
G}$. The best-fit model parameters of CoNb2O6 and $B\\!=\\!7$ T
($H\\!=\\!0.8H_{c}$) are used.

Generally, a magnon with the momentum ${\bf k}$ is kinematically allowed to
decay into a pair of magnons if the energy conservation for the single-magnon
energy and that of the two-magnon continuum

$\displaystyle\bar{\varepsilon}_{\bf k}=E_{{\bf k},{\bf q}},$ (87)

can be satisfied for some momenta ${\bf q}$. This condition naturally
separates the ${\bf k}$-space into the stable one, in which (87) cannot be
fulfilled, and the decay region, in which it is fulfilled [46]. They are easy
to visualize as having or not having an overlap of the single-magnon branch
with the two-magnon continuum, see Fig. 17(a), with the decay threshold
boundaries ${\bf k}^{*}$ and ${\bf G}-{\bf k}^{*}$ separating the stable
region from the decay region.

Such a threshold, or the entry-point of the single-magnon branch into the two-
magnon continuum, is necessarily a crossing of the single-magnon branch with
the minimum of the two-magnon continuum at that ${\bf k}^{*}$, which, in turn,
must correspond to a singularity in it as is discussed above, see also Refs.
[71, 46].

In the case relevant to the 1D model of CoNb2O6 discussed in this work, the
singularity associated with the minimum of the two-magnon continuum away from
the ${\bf\Gamma}$ point corresponds to the equivalent-magnon solution of the
extremum condition in (85), yielding an equation on the decay threshold
boundary ${\bf k}^{*}$

$\displaystyle\bar{\varepsilon}_{{\bf k}^{*}}=2\bar{\varepsilon}_{({\bf
k}^{*}-{\bf Q})/2},$ (88)

for $\Gamma\\!\lesssim\\!{\bf k}^{*}\\!\lesssim\\!\Gamma^{\prime}$. While this
equation can be solved numerically for the actual magnon energy
$\bar{\varepsilon}_{\bf k}$ of the considered model, a simplified solution for
${\bf k}^{*}$ can be derived analytically using the nearest-neighbor hopping
approximation $\bar{\varepsilon}_{\bf k}\\!\approx\\!{\sf E_{0}}+{\sf
J_{1}}\gamma_{\bf k}^{(1)}$ discussed above, which closely describes
$\bar{\varepsilon}_{\bf k}$ when the gap $\Delta_{0}$ is not too small,
yielding

$\displaystyle{\bf k}^{*}=2\arcsin\left(\frac{\sqrt{5+2\Delta}-1}{2}\right),$
(89)

where $\Delta\\!=\\!\Delta_{0}/|{\sf J_{1}}|\\!=\\!{\sf E_{0}}/|{\sf
J_{1}}|-1$ is the dimensionless gap.

One should also add that the $1/\sqrt{|\Delta{\bf k}|}$ singularities in the
on-shell magnon spectrum discussed in Sec. IV.2 are naturally connected to the
$1/\sqrt{|\Delta\omega|}$ one-dimensional Van Hove singularities in the two-
magnon DoS, with the latter transferred onto the single-magnon energies via
the anharmonic coupling.

The kinematic consideration of the energy conservation in Eq. (87) also allows
to trace the evolution of the decay surfaces, that is, the sets of the ${\bf
q}$ values into which decays are possible, as one traverses along the
$\bar{\varepsilon}_{\bf k}$ curve vs ${\bf k}$. In the 1D problem considered
here, the decay surfaces are, in fact, the pairs of the discrete ${\bf
q}$-points. At the decay threshold ${\bf k}^{*}$, they correspond to the two
equivalent points ${\bf q}^{*}$ and ${\bf q}^{*}+{\bf G}$, as one can see in
Fig. 17(b). As a function of ${\bf k}$, this set of pairs of the ${\bf
q}$-points traces a continuous line of the elliptic shape shown in Fig. 17(b),
with the schematics in Fig. 17(a) illustrating the decay process. The best-fit
parameters of CoNb2O6 and $B\\!=\\!7$ T ($H\\!=\\!0.8H_{c}$) are used in both
figures.