# Optimization of RIS Configurations for Multiple-RIS-Aided mmWave Positioning
Systems based on CRLB Analysis

Yu Liu, Sheng Hong, Cunhua Pan, Yinlu Wang, Yijin Pan, and Ming Chen
(Corresponding author: Cunhua Pan)C. Pan is with the School of Electronic
Engineering and Computer Science at Queen Mary University of London, London E1
4NS, U.K. (e-mail: c.pan@qmul.ac.uk).Y. Liu, Y. Wang, Y. Pan and M. Chen are
with the National Mobile Communications Research Laboratory, Southeast
University, Nanjing 210096, China. (e-mail:{liuyu_1994, yinluwang, chenming,
panyj}@seu.edu.cn).S. Hong is with School of Information Engineering at
Nanchang University, Nanchang 330100, China. (e-mail<EMAIL_ADDRESS>
###### Abstract

Reconfigurable intelligent surface (RIS) is a promising technology for future
millimeter-wave (mmWave) communication systems. However, its potential
benefits of adopting RIS for high-precision positioning in mmWave systems are
still less understood. In this paper, we study a multiple-RIS-aided mmWave
positioning system and derive the Cram$\rm{\acute{e}}$r-Rao error bound. Based
on the derived bound, we optimize the phase shift of the RISs by the particle
swarm optimization (PSO) algorithm. Numerical results have demonstrated the
advantages of using multiple RISs in enhancing the positioning accuracy in
mmWave systems.

###### Index Terms:

RIS, IRS, mmWave, MIMO, CRB, positioning

##  I Introduction

In the current 4G/5G wireless networks, local positioning system (LPS) has
attracted extensive research interests due to its key role in providing
accurate beamforming [1]. Recent research efforts in LPS have successfully
reduced the positioning error to the centimeter level with a single base
station in mmWave communications [2]. However, this level of accuracy cannot
meet the requirements demanded by emerging applications in future B5G/6G
networks, such as smart factories, autonomous driving, and augmented reality.
In these scenarios, the wireless propagation links suffer from severe
propagation loss and obstacles blockage when operating in mmWave and Terahertz
frequency band, which will deteriorate the positioning accuracy and
reliability [3, 4]. Therefore, it is imperative to improve the positioning
accuracy and reliability of LPS for the next-generation wireless networks.

Fortunately, the promising reconfigurable intelligent surface (RIS) technique
could provide a stable and high-precision positioning performance in a cost-
effective and energy-efficient manner. An RIS is a planar surface composed of
a large number of low-cost and passive reflecting elements, which can adjust
the amplitudes and phase shifts of the incident signals [5, 6, 7, 8]. Hence,
RIS can be deployed to create an alternative link to bypass the blockage of
the direct link, which enables the LPS to recover positioning functionality in
dead zones. Unlike fixed scatters in the environment, the RIS could be
envisioned as a controllable scatter by adjusting the phase shifts of the
reflection elements. Furthermore, an RIS with large aperture provides a higher
beamforming gain as well as a higher angle resolution, which is appealing for
positioning. It is expected that RIS could improve the positioning accuracy,
which is measured by the Cram$\rm{\acute{e}}$r-Rao lower bound (CRLB).

The authors in [9] investigated the advantages of deploying a single RIS with
a uniform linear array (ULA) for a two-dimensional positioning system. This
work was extended to a three-dimensional scenario with a planar RIS in [10].
However, to the best of our knowledge, there are no existing works that
adopted multiple RISs with varying sizes for positioning in mmWave systems. In
addition, the existing works designed the phase shifts of the RIS using the
beam matching method, which is a suboptimal choice in terms of minimizing the
CRLB.

The contributions of this paper are summarized as follows:

  1. 1.

We derive the position and orientation error bounds from the CRLB of a
multiple-RIS-aided positioning system.

  2. 2.

Due to the complicated expression of CRLB, we implement the heuristic
algorithm like particles swarm optimization (PSO) algorithm to achieve the
nearly optimal solution.

  3. 3.

The impact of various system parameters, such as the number of the RISs and
the size of the RIS, on the localization error bounds is studied in the
simulations.

##  II System Model

We consider a multiple-RIS-aided 3D positioning system shown in Fig. 1, which
consists of a base station (BS) equipped with a ULA of $N_{t}$ antennas, a
mobile user (MU) with a ULA of $N_{r}$ antennas and $K$ RISs with a uniform
planar array (UPA) of $L^{2}$ reflecting elements. The locations of the BS
antenna array, the $k$-th RIS, and the MU are respectively denoted by
$\mathbf{q}=\left[q_{x},q_{y},q_{z}\right]^{T}\in R^{3}$,
$\mathbf{s}_{k}=\left[s_{kx},s_{ky},s_{kz}\right]^{T}\in R^{3}$ and
$\mathbf{p}=\left[p_{x},p_{y},0\right]^{T}\in R^{3}$, where $q_{k}$ and
$s_{kz}$ symbolizes the heights of BS and the $k$-th RIS relative to the MU on
the ground. The rotation angle of the MU’s antenna array is denoted by
$\alpha\in\left(0,\pi\right]$, and no rotation is assumed for BS and RIS
arrays. The locations of the BS and the RISs are assumed to be known. Thus,
the objective of the positioning system is to estimate the coordinates and
rotation angle of the MU.

Figure 1: The positioning simulation scenario with three RISs

###  II-A Transmitter Model

We consider a system operating at a carrier frequency $f_{c}$ with wavelength
$\lambda$. The total bandwidth is $B$ with $N$ orthogonal frequency division
multiplexing (OFDM) subcarriers. The transmit signal on the $n$-th subcarrier
is denoted by $\mathbf{x}[n]=[x_{1}[n],\cdots,x_{M_{t}}[n]]^{T}$, where
$M_{t}\ll N_{t}$ denotes the number of beams. To simplify the analysis, it is
assumed that the pilot signal has unit power and flat spectrum. Thus, its
power spectrum is a constant of $\left|X(\omega)\right|^{2}=\frac{1}{2\pi B}$
within $\omega\in\left[0,2\pi B\right]$. The transmitted signals over
subcarrier $n$ can be expressed as $\bm{F}\mathbf{x}[n]$, where $\mathbf{F}$
is the beamforming matrix defined as
$\mathbf{F}=\left[\mathbf{f}_{1},\mathbf{f}_{2},\cdots,\mathbf{f}_{M_{t}}\right]\in
C^{N_{t}\times M_{t}}$, with $\mathbf{f}_{i}\in C^{N_{t}}$ denoting the unit-
norm transmitting vector.

###  II-B Wireless Channel Model

The channel state information (CSI) matrix $\mathbf{H}$ between the BS and the
MU could be expressed as the combination of the line-of-sight (LoS) component
and the reflecting component, which is given by

$\bm{H}[n]=\underbrace{\bm{H}_{0}[n]}_{\text{LoS
channel}}+\underbrace{\sum_{k=1}^{K}\bm{H}_{k}[n]}_{\text{reflecting
channel}}.$ (1)

The expressions of $\bm{H}_{0}[n]$ and $\bm{H}_{k}[n]$ are respectively given
by

$\displaystyle\bm{H}_{0}[n]\\!=\\!\gamma_{0}h_{0}\bm{\alpha}_{RX}\left(\theta_{RX,0}\right)\bm{\alpha}_{TX}^{H}\left(\theta_{TX,0}\right)e^{j2\pi
B\frac{n}{N}\tau_{0}},$ (2)
$\displaystyle\bm{H}_{k}[n]\\!=\\!\gamma_{k}h_{k}\bm{H}_{IM,k}\bm{\Theta}_{k}\bm{H}_{BI,k}e^{j2\pi
B\frac{n}{N}\tau_{k}},$ (3)

where $\gamma_{k}=\sqrt{N_{t}N_{r}/\rho_{k}}$ with $\rho_{k}$ being the path
loss of the $k$-th reflecting channel, $\tau_{k}$ is the propagation delay,
and $h_{k}$ is the complex channel gain. In Eq. (3), $\bm{H}_{BI,k}$ is the
channel matrix from the BS to the $k$-th RIS, and $\bm{H}_{IM,k}$ is the
channel matrix from the $k$-th RIS to the MU, and $\bm{\Theta}_{k}$ is the
diagonal phase shift matrix of the $k$-th RIS. They are respectively given by

$\displaystyle\bm{H}_{BI,k}=\bm{\alpha}_{RIS,IN}\left(\varphi_{in,k}^{a},\varphi_{in,k}^{e}\right)\bm{\alpha}_{TX}^{H}\left(\theta_{TX,k}\right),$
(4)
$\displaystyle\bm{H}_{IM,k}=\bm{\alpha}_{RX}\left(\theta_{RX,k}\right)\bm{\alpha}_{RIS,OUT}^{H}\left(\varphi_{out,k}^{a},\varphi_{out,k}^{e}\right),$
(5) $\displaystyle\bm{\Theta}_{k}=\delta\times
diag(e^{j\theta_{1}},e^{j\theta_{2}},\cdots,e^{j\theta_{L^{2}}})\in
C^{L^{2}\times L^{2}}.$ (6)

In Eq. (4) and Eq. (5),
$\bm{\alpha}_{TX}\left(\theta_{TX,k}\right)\in\mathbb{C}^{N_{t}}$ and
$\bm{\alpha}_{RX}\left(\theta_{RX,k}\right)\in\mathbb{C}^{N_{r}}$ are the
antenna response vectors of the transmitter and the receiver respectively,
with $\theta_{TX,k}$ denoting the angle of departure (AOD) and $\theta_{RX,k}$
the angle of arrival (AOA). The $i$-th entry of
$\bm{\alpha}_{TX}\left(\theta_{TX}\right)$ is
$\left[\bm{\alpha}_{TX}\left(\theta_{TX}\right)\right]_{i}=e^{j(i-1)\frac{2\pi}{\lambda}d\sin(\theta_{TX})}$
and the $i$-th entry of $\bm{\alpha}_{RX}\left(\theta_{RX}\right)$ is
$\left[\bm{\alpha}_{RX}\left(\theta_{RX}\right)\right]_{i}=e^{j(i-1)\frac{2\pi}{\lambda}d\sin(\theta_{RX})}$
, where $d$ is the distance between adjacent antennas and $\lambda=c/f_{c}$
with $c$ being the speed of light. Also in Eq. (4) and Eq. (5),
$\bm{\alpha}_{RIS,IN}(\varphi_{in,k}^{a},\varphi_{in,k}^{e})$ and
$\bm{\alpha}_{RIS,OUT}(\varphi_{out,k}^{a},\varphi_{out,k}^{e})$ are the array
response vectors of the BS-RISk and RISk-MU links respectively, where
$\varphi_{in,k}^{a}$ and $\varphi_{in,k}^{e}$ are respectively the azimuth AOA
and elevation AOA at the BS-RISk link, $\varphi_{out,k}^{a}$ and
$\varphi_{out,k}^{e}$ are respectively the azimuth AOD and elevation AOD at
the RISk-MU link. The $\left[a+(b-1)L\right]$-th elements in
$\bm{\alpha}_{RIS,IN}$ and $\bm{\alpha}_{RIS,OUT}$ are respectively given by

$\displaystyle\begin{aligned}
\left[\bm{\alpha}_{RIS,IN}\right]_{a+(b-1)L}=&e^{j\frac{2\pi}{\lambda}(b-1)d\cos\left(\varphi_{in}^{e}\right)}\\\
&\times
e^{j\frac{2\pi}{\lambda}d(a-1)\sin\left(\varphi_{in}^{e}\right)\sin\left(\varphi_{in}^{a}\right)},\\\
\end{aligned}$ (7) $\displaystyle\begin{aligned}
\left[\bm{\alpha}_{RIS,OUT}\right]_{a+(b-1)L}=&e^{j\frac{2\pi}{\lambda}(b-1)d\cos\left(\varphi_{out}^{e}\right)}\\\
&\times
e^{j\frac{2\pi}{\lambda}d(a-1)\sin\left(\varphi_{out}^{e}\right)\sin\left(\varphi_{out}^{a}\right)}.\end{aligned}$
(8)

In Eq. (6), $\delta=1$ represents the amplitudes of reflecting coefficients
and $\theta_{i}$ is the phase shift of the $i$-th reflecting element. Based on
[11], the path loss $\rho_{0}$ of the LoS channel is modeled as

$\displaystyle PL_{LOS}\left(d_{0}\right)=10\log_{10}\rho_{0}$ (9)
$\displaystyle=10\log_{10}(64\pi^{3})+10\alpha_{0}\log_{10}d_{0}+20\log_{10}f_{c}+\xi_{0},$

where $d_{0}$ is the distance between the MU and the BS, and $\alpha_{0}$ is
the path loss exponent,
$\xi_{0}\sim\mathcal{CN}\left(0,\sigma_{SF,0}^{2}\right)$ is the log-normal
shadow fading, with $\sigma_{SF,0}^{2}$ denoting the shadowing variance. The
path loss $\rho_{k}$ of the $k$-th RIS-aided link is modeled as

$\displaystyle PL_{RIS}\left(d_{k1},d_{k2}\right)=10\log_{10}\rho_{k}$ (10)
$\displaystyle=10\log_{10}(64\pi^{3})+10\alpha_{k}\log_{10}(d_{k1}\cdot
d_{k2})+40\log_{10}f_{c}+\xi_{k},$

where $d_{k1}$ is the distance of the BS-RIS link, $d_{k2}$ is the distance of
the RIS-MU link, $\alpha_{1}$ is the path loss exponent of the $k$-th
reflecting link, and $\xi_{k}\sim\mathcal{CN}\left(0,\sigma_{SF,k}^{2}\right)$
is the log-normal shadow fading, with $\sigma_{SF,k}^{2}$ denoting the
shadowing variance.

###  II-C Receiver Model

The received signal at the MU is given by

$\bm{y}\left[n\right]=\sqrt{P_{TX}}\bm{H}[n]\bm{F{\rm x}}[n]+\bm{{\rm n}}[n],$
(11)

where $P_{TX}$ is the transmitter power and $\bm{{\rm
n}}\sim\mathcal{N}(0,N_{0}\bm{I}_{N_{r}})$ is the additive white Gaussian
noise. In the following, we first derive the CRLB of the position and
orientation of the MU based on the received signal model in Eq. (11). Then, we
provide a low-complexity algorithm to optimize the phase shift of the RIS to
minimize the CRLB.

##  III Problem Formulation

We employ the two-stage approach to estimate the position and orientation of
the MU. In the first step, we estimate the channel parameters $\bm{\eta}$ in a
general form as

$\hat{\bm{\eta}}=\bm{\eta}+\bm{w},$ (12)

where $\bm{w}\sim\mathcal{CN}(\bm{H}[n]\bm{F{\rm x}}[n],\bm{\Sigma})$ denotes
the estimation error with $\bm{\Sigma}=N_{0}\bm{I}_{N_{r}}$ and $\bm{\eta}$
denotes the channel parameters, such as TOA, AOA, AOD and complex channel
gain. Specifically, $\bm{\eta}$ is given by

$\bm{\eta}=\left[\bm{\tau},\theta_{TX,0},\bm{\theta}_{RX},\bm{\varphi}_{out}^{a},\bm{\varphi}_{out}^{e},\bm{h}_{R},\bm{h}_{I}\right]^{T},$
(13)

where

$\displaystyle\bm{\tau}=\left[\tau_{0},\tau_{1},\cdots,\tau_{K}\right]^{T},$
(14)
$\displaystyle\bm{\theta}_{RX}=\left[\theta_{RX,0},\theta_{RX,1},\cdots,\theta_{RX,K}\right]^{T},$
(15)
$\displaystyle\bm{\varphi}_{out}^{a}=\left[\varphi_{out,1}^{a},\varphi_{out,2}^{a},\cdots,\varphi_{out,K}^{a}\right]^{T},$
(16)
$\displaystyle\bm{\varphi}_{out}^{e}=\left[\varphi_{out,1}^{e},\varphi_{out,2}^{e},\cdots,\varphi_{out,K}^{e}\right]^{T},$
(17)
$\displaystyle\bm{h}_{R}=\left[h_{R,0},h_{R,1},\cdots,h_{R,K}\right]^{T},$
(18)
$\displaystyle\bm{h}_{I}=\left[h_{I,0},h_{I,1},\cdots,h_{I,K}\right]^{T}.$
(19)

In the second step, we estimate the MU’s position $\hat{\bm{p}}$ and
orientation $\hat{\alpha}$ as follows

$\displaystyle\left[\hat{\bm{p}},\hat{\alpha}\right]$
$\displaystyle=\mathop{\arg\max}_{\left[\bm{p},\alpha\right]}\quad
p\left(\hat{\bm{\eta}}|\bm{\eta}\left(\bm{p},\alpha\right)\right)$ (20)
$\displaystyle=\mathop{\arg\max}_{\left[\bm{p},\alpha\right]}\quad\left(\hat{\bm{\eta}}-\bm{\eta}\left(\bm{p},\alpha\right)\right)^{T}\bm{\Sigma}^{-1}\left(\hat{\bm{\eta}}-\bm{\eta}\left(\bm{p},\alpha\right)\right),$

where $\bm{\eta}\left(\bm{p},\alpha\right)$ is the function characterizing the
relationship between $\bm{\eta}$ and $\left[\bm{p},\alpha\right]$. In
practice, $\bm{\eta}$ could be obtained via some compressive sensing
techniques due to the sparsity of mmWave channels.

The goal of the multiple-RIS-aided positioning system is to minimize the
expected estimation error of the MU’s position and orientation, which is
measured as

$\displaystyle\rm{var}\left(\hat{\boldsymbol{p}}\right)=\mathbb{E}\left[\left(\boldsymbol{p}-\hat{\boldsymbol{p}}\right)^{T}\left(\boldsymbol{p}-\hat{\boldsymbol{p}}\right)\right],$
(21)
$\displaystyle\rm{var}\left(\alpha\right)=\mathbb{E}\left[\left(\alpha-\hat{\alpha}\right)^{2}\right].$
(22)

##  IV Derivation of the Fundamental Bounds

In this section, we derive the CRLBs of the estimation of the MU’s position
and the rotation angle, based on which we can obtain the position error bound
(PEB) and rotation error bound (REB). To this end, we first derive the fisher
information matrix (FIM) from the channel parameters $\bm{\eta}$. The FIM of
the aforementioned parameters is a $(6K+5)\times(6K+5)$ matrix $\bm{J_{\eta}}$
with the following expression

$\left[\bm{J}_{\bm{\eta}}\right]_{mn}=\frac{2P_{TX}}{N_{0}B}\sum_{n=1}^{N}\mathrm{Re}\left\\{\frac{\partial\bm{\mu}^{H}[n]}{\partial\eta_{m}}\frac{\partial\bm{\mu}[n]}{\partial\eta_{n}}\right\\},$
(23)

where $\bm{\mu}[n]=\bm{H}[n]\bm{F{\rm x}}[n]$ since it is an OFDM system. The
components in $\bm{J}_{\bm{\eta}}$ are specified as follows:

$\displaystyle\frac{\partial\bm{\mu}[n]}{\partial\tau_{0}}=\gamma_{0}h_{0}\bm{\alpha}_{RX,0}\bm{\alpha}^{H}_{TX,0}\bm{F{\rm
x}}[n]e^{j2\pi B\frac{n}{N}\tau_{0}}j2\pi B\frac{n}{N},$ (24)
$\displaystyle\begin{aligned}
\frac{\partial\bm{\mu}[n]}{\partial\tau_{k}}=&\gamma_{k}h_{k}\bm{\alpha}_{RX,k}\bm{\alpha}_{OUT,k}^{H}\bm{\Theta}_{k}\bm{\alpha}_{IN,k}\bm{\alpha}^{H}_{TX,k}\\\
&\bm{F{\rm x}}[n]e^{j2\pi B\frac{n}{N}\tau_{k}}j2\pi
B\frac{n}{N},\end{aligned}$ (25)
$\displaystyle\frac{\partial\bm{\mu}[n]}{\partial\theta_{TX,0}}=\gamma_{0}h_{0}\bm{\alpha}_{RX,0}\bm{\alpha}^{H}_{TX,0}\bm{D}^{H}_{TX,0}\bm{F{\rm
x}}[n],$ (26)
$\displaystyle\frac{\partial\bm{\mu}[n]}{\partial\theta_{RX,0}}=\gamma_{0}h_{0}\bm{D}_{RX,0}\bm{\alpha}_{RX,0}\bm{\alpha}^{H}_{TX,0}\bm{F{\rm
x}}[n],$ (27) $\displaystyle\begin{aligned}
\frac{\partial\bm{\mu}[n]}{\partial\theta_{RX,k}}=&\gamma_{k}h_{k}\bm{D}_{RX,k}\bm{\alpha}_{RX,k}\bm{\alpha}_{OUT,k}^{H}\bm{\Theta}_{k}\\\
&\bm{\alpha}_{IN,k}\bm{\alpha}^{H}_{TX,k}\bm{F{\rm x}}[n],\\\ \end{aligned}$
(28) $\displaystyle\begin{aligned}
\frac{\partial\bm{\mu}[n]}{\partial\varphi_{out,k}^{a}}=&\gamma_{k}h_{k}\bm{\alpha}_{RX,k}\bm{\alpha}_{OUT,k}^{H}{\rm
diag}\left(\bm{c}^{a}_{out,k}\right)^{H}\\\
&\bm{\Theta}_{k}\bm{\alpha}_{IN,k}\bm{\alpha}^{H}_{TX,k}\bm{F{\rm
x}}[n],\end{aligned}$ (29) $\displaystyle\begin{aligned}
\frac{\partial\bm{\mu}[n]}{\partial\varphi_{out,k}^{e}}=&\gamma_{k}h_{k}\bm{\alpha}_{RX,k}\bm{\alpha}_{OUT,k}^{H}{\rm
diag}\left(\bm{c}^{e}_{out,k}\right)^{H}\\\
&\bm{\Theta}_{k}\bm{\alpha}_{IN,k}\bm{\alpha}^{H}_{TX,k}\bm{F{\rm
x}}[n],\end{aligned}$ (30) $\displaystyle\frac{\partial\bm{\mu}[n]}{\partial
h_{R,0}}=\gamma_{0}\bm{\alpha}_{RX,0}\bm{\alpha}^{H}_{TX,0}\bm{F{\rm
x}}[n]e^{j2\pi B\frac{n}{N}\tau_{0}},$ (31)
$\displaystyle\frac{\partial\bm{\mu}[n]}{\partial
h_{I,0}}=j\gamma_{0}\bm{\alpha}_{RX,0}\bm{\alpha}^{H}_{TX,0}\bm{F{\rm
x}}[n]e^{j2\pi B\frac{n}{N}\tau_{0}},$ (32) $\displaystyle\begin{aligned}
\frac{\partial\bm{\mu}[n]}{\partial
h_{R,k}}=&\gamma_{k}\bm{\alpha}_{RX,k}\bm{\alpha}_{OUT,k}^{H}\bm{\Theta}_{k}\bm{\alpha}_{IN,k}\bm{\alpha}^{H}_{TX,k}\\\
&\bm{F{\rm x}}[n]e^{j2\pi B\frac{n}{N}\tau_{k}},\end{aligned}$ (33)
$\displaystyle\begin{aligned} \frac{\partial\bm{\mu}[n]}{\partial
h_{R,k}}=&j\gamma_{k}\bm{\alpha}_{RX,k}\bm{\alpha}_{OUT,k}^{H}\bm{\Theta}_{k}\bm{\alpha}_{IN,k}\bm{\alpha}^{H}_{TX,k}\\\
&\bm{F{\rm x}}[n]e^{j2\pi B\frac{n}{N}\tau_{k}},\end{aligned}$ (34)

where

$\displaystyle\bm{\alpha}_{TX,k}=\bm{\alpha}_{TX}\left(\theta_{TX,k}\right),\quad\bm{\alpha}_{RX,k}=\bm{\alpha}_{RX}\left(\theta_{TX,k}\right),$
(35)
$\displaystyle\bm{\alpha}_{IN,k}=\bm{\alpha}_{RIS,IN}\left(\varphi_{in,k}^{a},\varphi_{in,k}^{e}\right),$
$\displaystyle\bm{\alpha}_{OUT,k}=\bm{\alpha}_{RIS,OUT}\left(\varphi_{out,k}^{a},\varphi_{out,k}^{e}\right),$
$\displaystyle\\!\\!\\!\\!\\!\bm{D}_{TX,k}(\theta)=j\frac{2\pi}{\lambda}d\cos(\theta_{TX,k}){\rm
diag}\left(0,1,\cdots,N_{t}-1\right),$
$\displaystyle\\!\\!\\!\\!\\!\bm{D}_{RX,k}(\theta)=j\frac{2\pi}{\lambda}d\cos(\theta_{RX,k}){\rm
diag}\left(0,1,\cdots,N_{r}-1\right),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\left[\bm{c}_{out,k}^{a}\right]_{a+(b-1)L}=j\frac{2\pi}{\lambda}(a-1)d\cos(\varphi_{out,k}^{a})\sin(\varphi_{out,k}^{e}),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\begin{aligned}
\left[\bm{c}_{out,k}^{e}\right]_{a+(b-1)L}=&j\frac{2\pi}{\lambda}d[(a-1)\sin\left(\varphi_{out,k}^{a}\right)\\\
&\cos\left(\varphi_{out,k}^{e}\right)-(b-1)\sin\left(\varphi_{out,k}^{e}\right)].\end{aligned}$

Next, the bijective matrix $\bm{T}$ is given by

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\bm{T}\\!=\\!\begin{bmatrix}\frac{\partial\bm{\tau}}{\partial
p_{x}}&\frac{\partial\theta_{TX,0}}{\partial
p_{x}}&\frac{\partial\bm{\theta}_{RX}}{\partial
p_{x}}&\frac{\partial\bm{\varphi}_{out}^{a}}{\partial
p_{x}}&\frac{\partial\bm{\varphi}_{out}^{e}}{\partial
p_{x}}&\frac{\partial\bm{h}_{R}}{\partial
p_{x}}&\frac{\partial\bm{h}_{I}}{\partial p_{x}}\\\
\frac{\partial\bm{\tau}}{\partial p_{y}}&\frac{\partial\theta_{TX,0}}{\partial
p_{y}}&\frac{\partial\bm{\theta}_{RX}}{\partial
p_{y}}&\frac{\partial\bm{\varphi}_{out}^{a}}{\partial
p_{y}}&\frac{\partial\bm{\varphi}_{out}^{e}}{\partial
p_{y}}&\frac{\partial\bm{h}_{R}}{\partial
p_{y}}&\frac{\partial\bm{h}_{I}}{\partial p_{y}}\\\
\frac{\partial\bm{\tau}}{\partial\alpha}&\frac{\partial\theta_{TX,0}}{\partial\alpha}&\frac{\partial\bm{\theta}_{RX}}{\partial\alpha}&\frac{\partial\bm{\varphi}_{out}^{a}}{\partial\alpha}&\frac{\partial\bm{\varphi}_{out}^{e}}{\partial\alpha}&\frac{\partial\bm{h}_{R}}{\partial\alpha}&\frac{\partial\bm{h}_{I}}{\partial\alpha}\end{bmatrix}.$
(36)

By utilizing the following geometric relationships among the BS, RISs and the
MU shown in Fig. 1

$\displaystyle\tau_{0}=\left\|\bm{q}-\bm{p}\right\|/c,\quad\tau_{k}=\left\|\bm{q}-\bm{s}_{k}\right\|/c+\left\|\bm{p}-\bm{s}_{k}\right\|/c,$
(37)
$\displaystyle\theta_{TX,0}=\arcsin\left(\left(p_{x}-q_{x}\right)/\left\|\bm{q}-\bm{p}\right\|_{2}\right),$
(38)
$\displaystyle\varphi_{out,k}^{a}=\arcsin\left(\frac{p_{y}-s_{ky}}{\sqrt{\left(p_{x}-s_{kx}\right)^{2}+\left(p_{y}-s_{ky}\right)^{2}}}\right),$
(39)
$\displaystyle\varphi_{out,k}^{e}=\arccos\left(s_{kz}/\left\|\bm{p}-\bm{s}_{k}\right\|_{2}\right),$
(40)
$\displaystyle\theta_{RX,0}=\arcsin\left(\frac{\left(p_{x}-q_{x}\right)\cos\alpha-\left(p_{y}-q_{y}\right)\sin\alpha}{\left\|\bm{q}-\bm{p}\right\|_{2}}\right),$
(41)
$\displaystyle\\!\\!\\!\\!\\!\\!\theta_{RX,k}=\arcsin\left(\frac{\left(p_{x}-s_{kx}\right)\cos\alpha-\left(p_{y}-s_{ky}\right)\sin\alpha}{\left\|\bm{p}-\bm{s}_{k}\right\|_{2}}\right),$
(42)

the non-zero entries of the bijective transformation matrix $\bm{T}$ are given
by

$\displaystyle\frac{\partial\tau_{0}}{\partial
p_{x}}=\frac{1}{c}\frac{p_{x}-q_{x}}{\left\|\bm{q}-\bm{p}\right\|_{2}},\quad\frac{\partial\tau_{0}}{\partial
p_{y}}=\frac{1}{c}\frac{p_{y}-q_{y}}{\left\|\bm{q}-\bm{p}\right\|_{2}},$ (43)
$\displaystyle\frac{\partial\tau_{k}}{\partial
p_{x}}=\frac{1}{c}\frac{p_{x}-s_{kx}}{\left\|\bm{p}-\bm{s}_{k}\right\|_{2}},\quad\frac{\partial\tau_{k}}{\partial
p_{y}}=\frac{1}{c}\frac{p_{y}-s_{ky}}{\left\|\bm{p}-\bm{s}_{k}\right\|_{2}},$
(44) $\displaystyle\frac{\partial\theta_{TX,0}}{\partial
p_{x}}=\frac{\sqrt{\left\|p-q\right\|_{2}^{2}-\left(p_{x}-q_{x}\right)^{2}}}{\left\|\bm{q}-\bm{p}\right\|_{2}^{2}},$
(45) $\displaystyle\frac{\partial\theta_{TX,0}}{\partial
p_{y}}=\frac{\left(p_{x}-q_{x}\right)\left(p_{y}-q_{y}\right)}{\sqrt{\left\|\bm{q}-\bm{p}\right\|_{2}^{2}-\left(p_{x}-q_{x}\right)^{2}}\cdot\left\|\bm{q}-\bm{p}\right\|_{2}^{2}},$
(46) $\displaystyle\frac{\partial\varphi_{out,k}^{a}}{\partial
p_{x}}=-\frac{p_{y}-s_{ky}}{\left(p_{x}-s_{kx}\right)^{2}+\left(p_{y}-s_{ky}\right)^{2}},$
(47) $\displaystyle\frac{\partial\varphi_{out,k}^{a}}{\partial
p_{y}}=\frac{p_{x}-s_{kx}}{\left(p_{x}-s_{kx}\right)^{2}+\left(p_{y}-s_{ky}\right)^{2}},$
(48) $\displaystyle\frac{\partial\varphi_{out,k}^{e}}{\partial
p_{x}}=\frac{s_{kz}\left(p_{x}-s_{kx}\right)}{\left\|p-s_{k}\right\|_{2}^{2}\sqrt{\left\|p-s_{k}\right\|_{2}^{2}-\beta_{IRS}^{2}}},$
(49) $\displaystyle\frac{\partial\varphi_{out,k}^{e}}{\partial
p_{y}}=\frac{s_{kz}\left(p_{y}-s_{ky}\right)}{\left\|p-s_{k}\right\|_{2}^{2}\sqrt{\left\|p-s_{k}\right\|_{2}^{2}-\beta_{IRS}^{2}}},$
(50) $\displaystyle\frac{\partial\theta_{RX,0}}{\partial
p_{x}}\\!=\\!\frac{\cos\alpha-\frac{\left(p_{x}-q_{x}\right)\left[\left(p_{x}-q_{x}\right)\cos\alpha-\left(p_{y}-q_{y}\right)\sin\alpha\right]}{\left\|p-q\right\|_{2}^{2}}}{\sqrt{\left\|\bm{q}-\bm{p}\right\|_{2}^{2}-\left[\left(p_{x}-q_{x}\right)\cos\alpha-\left(p_{y}-q_{y}\right)\sin\alpha\right]^{2}}},$
(51) $\displaystyle\frac{\partial\theta_{RX,0}}{\partial
p_{y}}\\!=\\!\frac{\sin\alpha+\frac{\left(p_{y}-q_{y}\right)\left[\left(p_{x}-q_{x}\right)\cos\alpha-\left(p_{y}-q_{y}\right)\sin\alpha\right]}{\left\|p-q\right\|_{2}^{2}}}{\sqrt{\left\|\bm{q}-\bm{p}\right\|_{2}^{2}-\left[\left(p_{x}-q_{x}\right)\cos\alpha-\left(p_{y}-q_{y}\right)\sin\alpha\right]^{2}}},$
(52)
$\displaystyle\frac{\partial\theta_{RX,0}}{\partial\alpha}\\!=\\!\frac{-\left(p_{x}-q_{x}\right)\sin\alpha-\left(p_{y}-q_{y}\right)\cos\alpha}{\sqrt{\left\|\bm{q}-\bm{p}\right\|_{2}^{2}-\left[\left(p_{x}-q_{x}\right)\cos\alpha-\left(p_{y}-q_{y}\right)\sin\alpha\right]^{2}}},$
(53) $\displaystyle\\!\\!\\!\\!\\!\\!\frac{\partial\theta_{RX,k}}{\partial
p_{x}}\\!=\\!\frac{\cos\alpha-\frac{\left(p_{x}-s_{kx}\right)\left[\left(p_{x}-s_{kx}\right)\cos\alpha-\left(p_{y}-s_{ky}\right)\sin\alpha\right]}{\left\|p-s_{k}\right\|_{2}^{2}}}{\sqrt{\left\|\bm{p}-\bm{s}_{k}\right\|_{2}^{2}-\left[\left(p_{x}-s_{kx}\right)\cos\alpha-\left(p_{y}-s_{ky}\right)\sin\alpha\right]^{2}}},$
(54) $\displaystyle\\!\\!\\!\\!\\!\\!\frac{\partial\theta_{RX,k}}{\partial
p_{y}}\\!=\\!\frac{\sin\alpha+\frac{\left(p_{y}-s_{ky}\right)\left[\left(p_{x}-s_{kx}\right)\cos\alpha-\left(p_{y}-s_{ky}\right)\sin\alpha\right]}{\left\|p-s_{k}\right\|_{2}^{2}}}{\sqrt{\left\|\bm{p}-\bm{s}_{k}\right\|_{2}^{2}-\left[\left(p_{x}-s_{kx}\right)\cos\alpha-\left(p_{y}-s_{ky}\right)\sin\alpha\right]^{2}}},$
(55)
$\displaystyle\\!\\!\\!\\!\\!\\!\frac{\partial\theta_{RX,k}}{\partial\alpha}\\!=\\!\frac{-\left(p_{x}-s_{kx}\right)\sin\alpha-\left(p_{y}-s_{ky}\right)\cos\alpha}{\sqrt{\left\|\bm{p}-\bm{s}_{k}\right\|_{2}^{2}-\left[\left(p_{x}-s_{kx}\right)\cos\alpha-\left(p_{y}-s_{ky}\right)\sin\alpha\right]^{2}}}.$
(56)

Finally, we obtain the FIM $\bm{J}$ in position domain through $\bm{J}_{\eta}$
and $\bm{T}$,

$\bm{J}=\bm{TJ_{\eta}T}^{H}.$ (57)

Thus, the PEB is derived as the root square of the trace of the first $2\times
2$ sub-matrix of $\bm{J}$

${\rm PEB}=\sqrt{\mathrm{tr}\left(\bm{J}_{1:2,1:2}^{-1}\right)},$ (58)

and the REB is given by the root square of the third diagonal entry of
$\bm{J}$

${\rm REB}=\sqrt{\bm{J}_{3,3}^{-1}}.$ (59)

##  V Phase Shift Optimization

Based on the derived CRLB expression, it is observed that the phase shifts of
the RIS have a great impact on the CRLB. In this paper, we aim to optimize the
phase shifts of the RIS to minimize the sum PEB and REB, which is formally
formulated as follows:

$\displaystyle\min_{\bm{\Theta}_{k}}$ $\displaystyle\rm{PEB+REB}$ (60)
$\displaystyle\rm{s.t.}$
$\displaystyle\bm{\Theta_{k}}=\delta\times\rm{diag}\left(e^{j\theta_{1}},e^{j\theta_{2}},\cdots,e^{j\theta_{L^{2}}}\right),$
$\displaystyle\theta_{i}\in[0,2\pi),\quad i=1,2,\cdots,L^{2}.$

However, the objective function in Problem (60) is too complex and the
conventional gradient method is not applicable. To address this issue, the
heuristic method such as PSO is adopted, which only needs to evaluate the
objective function value in each iteration rather than calculating the fisrt-
order derivative of the original objective function that entails high
computational complexity.

##  VI Simulation Results

Consider a positioning system with one direct path and three reflecting paths
separately going through three RISs. The BS is located at
$\rm{(0\,m,0\,m,40\,m)}$ and MU is located at $\rm{(90\,m,30\,m,0\,m)}$. The
locations of the RISs are $\rm{(60\,m,45\,m,15\,m)}$,$\rm{(50\,m,50\,m,5\,m)}$
and $\rm{(40\,m,20\,m,10\,m)}$. The number of reflecting elements at the RIS
is set as $L^{2}=16^{2}=256$. The numbers of transmitter and receiver antennas
are $N_{t}=32$ and $N_{r}=8$ respectively. The carrier frequency is
$f_{c}=4.9$ GHz and the number of subcarriers is $N=128$. The noise power
spectrum density is $N_{0}=-174$ dBm/Hz, bandwidth is $B=20\,$MHz, the path
loss exponent of the direct channel is $\alpha_{0}=3.7$, the path loss
exponent of the $k$-th RIS is $\alpha_{k}=2.2$ and shadow fading parameters of
direct path and reflecting path are respectively $\sigma_{SF,0}=4$ and
$\sigma_{SF,k}=7$.

###  VI-A The Impact of the Number of RISs

In Fig. 2, we compare the performance of our optimized phase shifts of the RIS
with the case when the phase shifts are randomly set. It is observed from Fig.
2 that the optimized phase shifts can achieve much better performance than the
randomly generated phase shifts in terms of both PEB and REB. In addition,
Fig. 2 shows that PEB and REB gradually decreases with the number of RISs,
which demonstrates the advantages of using multiple RISs in positioning.

###  VI-B The Impact of the Phase Shifts of RIS

Concerning the optimization of the phase shifts, one method is to set the
phase shifts that align with the channel vectors as in [9, 10]. Specifically,
the phase shifts of the RIS are set as follows:

$\displaystyle\bm{\Theta}_{k}^{\prime}=\delta\times{\rm diag}$
$\displaystyle\left\\{\bm{\alpha}_{RIS,IN}^{*}\left(\varphi_{in,k}^{a},\varphi_{in,k}^{e}\right)\right.$
(61)
$\displaystyle\left.\odot\bm{\alpha}_{RIS,OUT}\left(\varphi_{out,k}^{a},\varphi_{out,k}^{e}\right)\right\\}.$

However, this approach only maximizes part of the elements in the FIM, and
cannot ensure that each entry of the FIM is maximized, especially for those
elements containing $\frac{\partial\bm{\mu}[n]}{\varphi_{out,k}^{a}}$ and
$\frac{\partial\bm{\mu}[n]}{\varphi_{out,k}^{e}}$. On the other hand, the PSO
algorithm aims to search the globally optimal solution. From Fig. 2, we find
that when the system has only one RIS, the performance of the PSO-optimized
phase shifts are close to that of the beam-aligned (BA) phase shifts.
Furthermore, in the case of multiple RISs, although the REBs of two schemes
are still close, the PSO-optimized phase shifts performs much better than the
beam-aligned phase shifts in terms of the PEB performance.

(a) PEB

(b) REB

Figure 2: The PEB and REB with different number of paths and different kinds
of phase shifts

###  VI-C The Impact of the Size of RIS

Fig. 3 illustrates the impact of the size of RIS on the positioning
performance. It is readily seen that more reflecting elements could contribute
to a higher resolution, which results in a higher quality of positioning. It
is interesting to find that the improvement is significant when the size is
small and will reduce as the size increases.

(a) PEB

(b) REB

Figure 3: The PEB and REB with different size of RISs

##  VII Conclusion

In this paper, the Cram$\rm{\acute{e}}$r-Rao bound of a multiple-RIS-aided
mmWave MIMO positioning system has been derived. Through the
Cram$\rm{\acute{e}}$r-Rao bound of location and rotation angle estimation,
termed as PEB and REB respectively, we have investigated the impact of the
number, the sizes, and the phase shifts of RISs on positioning performance.
Simulation results have shown the improvement in PEB and REB by adopting
multiple RISs, indicating the great potential of RISs to help achieve high
positioning accuracy.

## References

  * [1] C. Laoudias, A. Moreira, S. Kim, S. Lee, L. Wirola, and C. Fischione, “A survey of enabling technologies for network localization, tracking, and navigation,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 4, pp. 3607–3644, 2018.
  * [2] A. Shahmansoori, G. E. Garcia, G. Destino, G. Seco-Granados, and H. Wymeersch, “Position and orientation estimation through millimeter-wave mimo in 5g systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 1822–1835, 2017.
  * [3] M. R. Akdeniz, Y. Liu, M. K. Samimi, S. Sun, S. Rangan, T. S. Rappaport, and E. Erkip, “Millimeter wave channel modeling and cellular capacity evaluation,” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1164–1179, 2014.
  * [4] Z. Abu-Shaban, X. Zhou, T. Abhayapala, G. Seco-Granados, and H. Wymeersch, “Error bounds for uplink and downlink 3d localization in 5g millimeter wave systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 8, pp. 4939–4954, 2018.
  * [5] C. Pan, H. Ren, K. Wang, J. F. Kolb, M. Elkashlan, M. Chen, M. Di Renzo, Y. Hao, J. Wang, A. L. Swindlehurst, X. You, and L. Hanzo, “Reconfigurable intelligent surfaces for 6g systems: Principles, applications, and research directions,” _IEEE Commun. Mag._ , vol. 59, no. 6, pp. 14–20, 2021.
  * [6] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. de Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [7] A. Elzanaty, A. Guerra, F. Guidi, and M.-S. Alouini, “Reconfigurable intelligent surfaces for localization: Position and orientation error bounds,” _IEEE Trans. Signal Process._ , vol. 69, pp. 5386–5402, 2021.
  * [8] H. Zhang, H. Zhang, B. Di, K. Bian, Z. Han, and L. Song, “Metalocalization: Reconfigurable intelligent surface aided multi-user wireless indoor localization,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2021.
  * [9] J. He, H. Wymeersch, L. Kong, O. Silvén, and M. Juntti, “Large intelligent surface for positioning in millimeter wave mimo systems,” in _2020 IEEE 91st Vehicular Technology Conference (VTC2020-Spring)_ , 2020, pp. 1–5.
  * [10] Y. Liu, E. Liu, R. Wang, and Y. Geng, “Reconfigurable intelligent surface aided wireless localization,” in _ICC 2021 - IEEE International Conference on Communications_ , 2021.
  * [11] W. Tang, M. Z. Chen, X. Chen, J. Y. Dai, Y. Han, M. Di Renzo, Y. Zeng, S. Jin, Q. Cheng, and T. J. Cui, “Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 421–439, 2021.


#  Truncated string-state-space approach and its application to non-integrable
spin 1/2 Heisenberg chain

Jiahao Yang (gbsn杨家豪)  Tsung-Dao Lee Institute, Shanghai Jiao Tong University,
Shanghai 201210, China  Jianda Wu (gbsn吴建达<EMAIL_ADDRESS>Tsung-Dao Lee
Institute, Shanghai Jiao Tong University, Shanghai 201210, China  School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Shanghai Branch, Hefei National Laboratory, Shanghai 201315, China

###### Abstract

By circumventing the difficulty of obtaining exact string state solutions to
Bethe ansatz equations, we devise a truncated string-state-space approach for
investigating spin dynamics in non-integrable spin-1/2 Heisenberg chain
subjected to staggered field at various magnetization. The obtained dynamical
spectra reveal a series of elastic peaks at integer multiples of the ordering
wavevector $Q$, indicating the presence of multi-$Q$ Bethe string states
within the ground state. The spectrum exhibits a separation between different
string continua as the strength of the staggered field increases at low
magnetization,reflecting the confinement of the Bethe strings. The approach
has successfully provided a unified string-state-based picture for
understanding spin dynamics observed in various phases of quasi-1D
antiferromagnet $\rm YbAlO_{3}$ Yang  _et al._ (2023a), and also paves a way
for investigating low-dimensional non-integrable systems.

Introduction. One-dimensional quantum systems, characterized by the presence
of exact solutions and unique integrability, offers fascinating arena to study
quantum many-body physics. A significant subset of these systems can be
analytically solved through Bethe ansatz, with notable examples including the
one-dimensional (1D) spin-1/2 XXZ model Jimbo and Miwa (1995); Franchini
(2017); Yang  _et al._ (2023b); He  _et al._ (2017), the Gaudin-Yang model
Gaudin (1967); Yang (1967); Guan  _et al._ (2013), and the Lieb-Liniger model
Lieb and Liniger (1963); Lieb (1963). Although the Bethe ansatz has paved the
way for determining eigenstates and eigenenergies of those systems, it has
been a long challenge for calculating their form factors and thus dynamical
response, which is partly tackled recently Li  _et al._ (2023); Guan  _et al._
(2013); Guan and Andrei (2018); Kitanine  _et al._ (1999, 2000). Empowered by
the theoretical development, the existence of celebrated many-body
quasiparticles, such as spinons Caux and Hagemans (2006); Caux  _et al._
(2008a, b); Castillo (2020) and strings Takahashi (1971); Kohno (2009); Yang
_et al._ (2019), has been experimentally discovered within the quasi-1D
materials Yang  _et al._ (2023b); Lake  _et al._ (2013); Wang  _et al._
(2018); Bera  _et al._ (2020); Wang  _et al._ (2019). These progresses witness
the cooperative effort of both theorists and experimentalists to unveil the
intricate nature of these exotic phenomena.

Bearing real materials in mind, it becomes crucial to ask how robust
integrable physics is against nonintegrable perturbations that may partially
or fully break the conservation laws of the integrable systems. This has
inspired extensive research focused on nonintegrable models such as spin-1/2
ladders Steinigeweg  _et al._ (2016); Batchelor  _et al._ (2005), chains with
a staggered magnetic field Essler  _et al._ (1997); Starykh and Balents
(2014), frustrated spin chains Bonča  _et al._ (1994), and dimerized spin
chains Heidrich-Meisner  _et al._ (2003); Keselman  _et al._ (2020). Most of
them are performed by effective field theory Essler  _et al._ (1997); Starykh
and Balents (2014) or the numerical methods such as exact diagonalization (ED)
Bonča  _et al._ (1994); Heidrich-Meisner  _et al._ (2003), matrix-product-
state Keselman  _et al._ (2020), and Quantum Monte Carlo methods Zhou  _et
al._ (2021). However, on one hand the numerical method in general lacks a
clear understanding of essential physical picture, on the other hand the
effective field theory can only provide a limited insight within the low-
energy and long-wavelength limit. Therefore, a method able to go beyond those
limitations is always desired.

At first glance it may be promising to apply Bethe states to study non-
integrable systems. However, a notorious open problem persists: finding the
precise complex solutions of the Bethe ansatz equation (BAE) for string states
Takahashi (1971); Takahashi and Suzuki (1972); Fujita  _et al._ (2003);
Ilakovac  _et al._ (1999); Isler and Paranjape (1993); Essler  _et al._
(1992). In the past decades, many approaches have been explored, including
carefully-designed iterative method Hagemans and Caux (2007) and rational
$Q$-system method Marboe and Volin (2017); Bajnok  _et al._ (2020); Hou  _et
al._ (2023). The former can easily access large system size but suffer from
many unphysical solutions with repeated roots. Although the later can solve
the BAE for all exact solutions simultaneously, it is limited to small system
size. Those shortcomings impede the practical application of the Bethe states
to understand non-integrable system with reasonable size.

In this manuscript, we first devise an algorithm for spin-1/2 Heisenberg
chain, which is able to obtain the exact solutions of Bethe string states.
Based on these states with proper eigenenergy cutoff, we develop a truncated
string-state-space approach (TS3A) to study the non-integrable spin dynamics
in spin-1/2 Heisenberg chain with staggered field. The TS3A can determine the
asymptotically exact eigenstates for non-integrable Hamiltonians, whose
accuracy can be systematically improved by raising the energy cutoff and the
number of types of string states.

Following the TS3A, we analyze the dynamic spectra of the non-integrable spin
dynamics in spin-1/2 Heisenberg chain with staggered field characterized by
wavevector $Q$. In addition to the $Q$-ordering of the system, a series of
elastic peaks appear at $nQ\ (n=2,3,\ldots)$, since the ground state contains
multi-$Q$ Bethe string states. Furthermore, the staggered field plays a role
of the confining field for the Bethe string states, leading to the confinement
effect in inelastic spectra. At low magnetization, the spin dynamics are
primarily governed by length-1 Bethe strings with small staggered field.
However, as staggered field increases, the dominance within the high-energy
sector of the spectra transfers to confined length-2 Bethe strings. We further
benchmark our results with ED calculation, and the comparison yields excellent
agreement. Notably, the TS3A has successfully provided a unified Bethe-string-
based understanding for the experimental observations in the quasi-1D
antiferromagnet $\rm YbAlO_{3}$ Yang  _et al._ (2023a).

Model. Our parent Hamiltonian is the 1D Heisenberg spin-1/2 model with
longitudinal field $h_{z}$,

$H_{0}=J\sum_{n=1}^{N}\mathbf{S}_{n}\cdot\mathbf{S}_{n+1}-h_{z}S_{n}^{z}$ (1)

with total sites $N$, antiferromagnetic coupling $J$, and spin operators
$\mathbf{S}_{n}$ with components $S_{n}^{\mu}$ ($\mu=x,y,z$) at site $n$. By
introducing a staggered field
$\mathbf{h}_{Q}=h_{Q}\sum_{i}\cos(Qr_{i})\hat{z}$ which couples to the spin
chain, the total Hamiltonian becomes non-integrable,

$H=H_{0}+H^{\prime},$ (2)

where

$H^{\prime}=-\sum_{n}\mathbf{h}_{Q}\cdot\mathbf{S}_{n}=-h_{Q}\sum_{n}\cos(Qr_{n})S^{z}_{n}.$
(3)

The $h_{Q}$ is the strength of the staggered field with ordering wavevector
$Q=(1-m)\pi$, where $m$ is the magnetization density, the ratio of
magnetization $M_{z}$ to its saturation value $M_{s}$. In practice, the
staggered field can be effectively induced from 3D magnetic ordering of
quasi-1D materials, such as $\rm YbAlO_{3}$ Wu  _et al._ (2019); Fan  _et al._
(2020); Nikitin  _et al._ (2021); Yang  _et al._ (2023a), $\rm
SrCo_{2}V_{2}O_{8}$ Bera  _et al._ (2014) and $\rm BaCo_{2}V_{2}O_{8}$ Zou
_et al._ (2021). We note that the staggered field can be both commensurate and
incommensurate cases, depending on whether $2\pi/Q$ is a rational or
irrational number, respectively.

For investigating spin dynamics of the system we focus on the zero-temperature
dynamical structure factor (DSF) for spin along longitudinal ($z$) direction
($\hbar=1$),

$D^{zz}(q,\omega)=2\pi\sum_{\mu}|\langle
GS|S_{q}^{z}|\mu\rangle|^{2}\delta(\omega-E_{\mu}+E_{GS}).$ (4)

with transfer momentum $q$ and transfer energy $\omega$ between the ground
state $|GS\rangle$ and excited states $|\mu\rangle$ whose energies are
$E_{GS}$ and $E_{\mu}$, respectively. The form factor $\langle
GS|S_{q}^{z}|\mu\rangle$ is accessible for the integrable Hamiltonian $H_{0}$
Kohno (2009); Caux and Hagemans (2006); Caux  _et al._ (2008b, a); Castillo
(2020); Yang  _et al._ (2019). In order to tackle non-integrable Hamiltonian
$H$ [Eq. (2)], we propose the TS3A to obtain asymptotically exact result for
DSF, as discussed in the following.

Figure 1:  (a) The rapidities of Bethe strings with different lengths in the
complex plane. (b-d) Pictorial spin configurations of the Bethe strings.
Viewing the down spin (red arrow) and up spin (blue arrow) as the magnon and
vacuum, respectively. A single Bethe string $\chi^{j}$ contains $j$ bounded
magnons and lives in the corresponding $j$-string band in (e). (e) The energy-
momentum relation for Bethe strings. (f) Illustration of truncated string-
state-space.

Bethe string state. To construct the Hilbert space of Bethe string states for
the TS3A , we begin with an introduction on the coordinate Bethe ansatz, and
an exposure of an efficient algorithm for obtaining the exact solutions from
the BAE.

Due to U(1) symmetry of $H_{0}$ [Eq. (1)] the magnetization $M_{z}=1/2-M/N$ is
the conserved quantity, where $M$ is the number of flipped spins with respect
to the fully up-polarized state. The flipped spins are viewed as particles,
viz, magnons. In the coordinate Bethe ansatz Bethe (1931); Karabach  _et al._
(1997); Karbach  _et al._ (2000), the eigenstate of $H_{0}$ [Eq. (1)] is the
Bethe state with $M$ magnons, which is determined by a set of rapidities
$\\{\lambda_{j}\\}_{M}$ satisfying the BAE,

$\left(\frac{\lambda_{j}+i/2}{\lambda_{j}-i/2}\right)^{N}=\prod_{k\neq
j}\left(\frac{\lambda_{j}-\lambda_{k}+i}{\lambda_{j}-\lambda_{k}-i}\right).$
(5)

These rapidities $\\{\lambda_{j}\\}_{M}$ manifest as either complex-conjugate
pairs or real numbers [Fig. 1(a)] Vladimirov (1986). The pair of complex
rapidities implies significant physical property: the corresponding magnons
exhibit an intriguing phenomenon in coordinate space, forming an effective
many-body bound magnons commonly referred to as a “Bethe string” Takahashi
(1971); Karabach  _et al._ (1997); Takahashi (1999). And the length of the
string is determined by the number of rapidities with a common real center.
Intuitively, Bethe string $\chi^{j}$ ($j\geq 2$) of length $j$ is a “big”
quasiparticle wherein $j$ bounded magnons move coherently, referred to as
$j$-string [Fig. 1(b-d)]. When $j=1$, the 1-string $\chi^{1}$ is just the
unbound magnon. Correspondingly, the rapidities of a string $\chi^{j}$ takes
the form Takahashi (1971); Takahashi and Suzuki (1972)

$\lambda^{n}_{j,\alpha}=\lambda_{j,\alpha}+\frac{i}{2}(j+1-2n)+d^{n}_{j,\alpha},$
(6)

where $n=1,\ldots,j$ denotes the $j$th magnon in the $j$-string. The number of
$j$-strings is denoted as $M_{j}$, and $\alpha=1,\ldots,M_{j}$ label different
$j$-strings with the same length $j$. Thus we have $\sum_{j}jM_{j}=M$ for a
$M$-magnon Bethe state. We refer to $\lambda_{j,\alpha}$ as the string center
which gives the real part of the $j$-string if the deviation
$d^{n}_{j,\alpha}$ is omitted. Under the assumption $d^{n}_{j,\alpha}=0$, we
obtain eigenenergy of a Bethe string state,
$E=\sum_{j,\alpha}\varepsilon_{j,\alpha}$, with
$\varepsilon_{j,\alpha}={-2jJ}/{(4\lambda_{j,\alpha}^{2}+j^{2})}$. Therefore,
we can exhibit the energy-(quasi)momentum relation for different strings in
Fig. 1(e) with $k=\pi-2\arctan(2\lambda)$.

To ensure clarity in terminology, we refer to a Bethe state with all $M$
magnons being $\chi^{1}$ as the 1-string state. For $j>1$, we classify an
$n*j$-string state with $M_{j}=n$ and $M_{1}=M-n*j$. When $n=1$, it is simply
referred to as the $j$-string state. And this convention can be consistently
extended to cover other cases.

In practice, we observe that $d^{n}_{j,\alpha}=0$ may not be a suitable
assumption for our TS3A. The underlying issue lies in the fact that imprecise
solutions fail to provide an accurate matrix representation of $H^{\prime}$
[Eq. (3)], resulting in substantial errors in the eigenstates for $H$ [Eq.
(2)] after diagonalization. Finding exact solutions to BAE [Eq. (5)] remains a
long-standing and challenging problem. Efforts such as the iterative method
used in the Ref. Hagemans and Caux (2007) is unable to generate all exact
solutions when lattice sites $N\geq 12$. The limitations arise from the
generation of repeated real rapidities in certain string states which are
physically prohibited com . We find the key to solve the problem is that the
repeated real rapidities actually form a complex conjugate pair with minor
imaginary part (typically $\lesssim 1/N$). To implement the observation into
the algorithm, we divide and conquer. For instance, for a 3-string state in
Table 1, the typical rapidity pattern is that three of them share a common
real part up to a finite deviation $d_{j,\alpha}^{n}$, while the remaining
rapidities are all real. However, when we encounter the repeated real
rapidities, usually involving the 3-string center and one real rapidity, we
introduce a small imaginary part to create a complex conjugate pair. This pair
and other rapidities are then treated as a new initial guess for the BAE, from
which we are able to efficiently obtain the exact solution Deguchi and Giri
(2015).

| $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$
---|---|---|---|---|---
unphysical Hagemans and Caux (2007) | $0.4955+0.9622i$ | $0.4955-0.9622i$ | 0.4458 | 0.4458 | 0.1803
physical | $0.4918+0.9615i$ | $0.4918+0.9615i$ | $0.4448+0.0188i$ | $0.4448-0.0188i$ | 0.1807
Table 1:  Rapidities of Bethe ansatz equaitons for $N=12$ and $M=5$. The
solutions presented in the first row are obtained from the Ref. Hagemans and
Caux (2007) and that in second row are obtained from our method. Note that in
first row, $\lambda_{3}$ and $\lambda_{4}$ are coincided, representing an
unphysical outcome.  Figure 2:  The absolute value of matrix entry
$\langle\mbox{1-str}|S_{\pi}^{z}|B_{b}\rangle$ where $|B_{b}\rangle$ ranges
from 1-string to 2*2-string states.

Truncated string-state-space approach. In many physical problems, it is
usually sufficient for us to focus on the low-energy subspace instead of the
complete Hilbert space James  _et al._ (2018). Here, the TS3A [Fig. 1g] we
developed can yield asymptotically exact results, which can be systematically
refined to a satisfactory level. From Fig. 2, it is evident that the form
factor for two-string states generally diminishes rapidly as the difference in
string length increases. Therefore, we can safely truncate the string state
space into relevant subspace in terms of string lengths as well as a proper
energy cutoff.

When considering TS3A to $H_{0}$ [Eq. (2)], we pick up those short-string-
length states with a proper energy cutoff $E^{B}_{cut}$. To a good
approximation, it’s enough to consider Bethe states with various combination
of 1-, 2-, and 3-strings for the zero-temperature DSF Kohno (2009); Yang  _et
al._ (2019). Furthermore, when considering a non-integrable perturbation, such
as $H^{\prime}$ [Eq. (3)], the Bethe string state is no longer the eigenstate.
Thus, we apply the TS3A to obtain new ground state and excited states before
the calculation of any physical quantities like DSFs. To do so, the first step
is to construct the matrix representation of the non-integrable Hamiltonian
$H$ [Eq. (2)] within the truncated space, $H_{tr}=\delta_{ab}E^{B}_{a}+\langle
B_{a}|H^{\prime}|B_{b}\rangle$ with $E^{B}_{a}\leq E^{B}_{cut}$. The dimension
of $H_{tr}$ is typically much less than $2^{N}$ and is controlled by the
energy cutoff $E^{B}_{cut}$ and types of Bethe string states. After the
diagnolization of $H_{tr}$, we obtain the asymptotically exact eigenstates and
eigenenergies for $H$, which enable us to further study the physical
quantities such as zero-temperature DSF Eq. (4). Then we discuss the results
in the following.

Figure 3: Zero temperature DSF of XXX chain with lattice size $N=48$. The DSFs
with $h_{Q}=0J$, $0.2J$, and $0.4J$ (from left to right), and magnetization
density $m$=12.5%, 25%, 50%, 75% (from top to bottom).  Figure 4: Zero
temperature DSF of XXX chain with lattice size $N=16$. The DSFs with
$h_{Q}=0J$, $0.2J$, and $0.4J$ (from left to right), and magnetization density
$m$=12.5%, 25%, 50%, 75% (from top to bottom).  Figure 5: The static structure
factor with $h_{Q}=0J$, $0.2J$, and $0.4J$, and magnetization density
$m$=12.5%, 25%, 50%, 75%.

Results and discussion. With the staggered field perturbation term
$H^{\prime}$ [Eq. (3)], the Hamiltonian $H$ [Eq. (2)] is not longer
integrable. Following the TS3A introduced above, we obtain the approximately
exact eigenstates and eigenenergies. Due to the fact that
$\sum_{n}\cos(Qr_{n})S^{z}_{n}\propto(S^{z}_{Q}+S^{z}_{-Q})$, the non-
vanishing matrix element of $H^{\prime}$ only appears between Bethe states
with momentum difference $\Delta q=\pm Q$. As a result, the ground state
consists of Bethe states with momenta $nQ\ (n=1,2,\ldots)$ adding to original
ground state. For static structure factor $D^{zz}(q,\omega=0)$ [Fig. 5], a
series of staggered-field induced peaks appear at $nQ$ (n=1, 2, $\cdots$),
manifesting the presence of multi-$Q$ Bethe states in the ground state. These
results are consistent with the experimental observation of $\rm YbAlO_{3}$
and the multiple fermion scattering mechanism of the Tomonaga-Luttinger liquid
theory Nikitin  _et al._ (2021).

The staggered field also serves as a confining field for the Heisenberg spin
chains Essler  _et al._ (1997); Starykh and Balents (2014). When $h_{Q}=0$,
the spectra show gapless excitations at $q=(1\pm m)\pi$, and the 2-string
states are hard to be distinguished from the broad continuum of 1-string
states. When $h_{Q}>0$, the energy gap not only appears near the elastic line,
but also gradually develops between the continua of 1- and 2-string states as
$h_{Q}$ increases. The induced gaps reflect the energy cost for the excitation
of Bethe strings, which is known as confinement of Bethe strings and has been
observed in $\rm YbAlO_{3}$ Yang  _et al._ (2023a). At small magnetization
($m\lesssim 25\%$), 2-string states are confined and becomes distinguishable
from the confined 1-string continuum, however, the 2-string continuum fades
away into higher energy window as $m$ increases.

We further compare our results with the ED calculations in Fig. 4. This
comparison reveals a remarkable consistency in the overall shape of the
spectra. Moreover, our larger size calculation not only reduces finite size
effect but also renders the characteristics of spectra more transparent and
discernible. Crucially, our TS3A offers two key advantages compared to the ED
method: First, it has higher efficiency, facilitating much larger systems
($N\gtrsim 50$); second, it naturally provides a unified Bethe-string-based
physical picture for understanding the underlying physics.

Conclusion. We exploit an efficient algorithm to find exact solutions for the
Bethe string states from the BAE of the spin-1/2 Heisenberg spin chain. Based
on the exact solutions we further develop the TS3A which enables us to
determine eigenstates and eigenenergies of non-integrable spin-1/2 Heisenberg
systems with U(1) symmetry preserved. The method is then applied to
systematically study the dynamics of non-integrable spin-1/2 Heisenberg spin
chain under staggered field. In the dynamic spectra, we reveal a series of
elastic peaks locating at the integer multiples of the ordering wavevector
$Q$, signifying the existence of multi-$Q$ Bethe string states within the
ground state. Moreover, the staggered field serves as a confining field for
Bethe string states, giving rise to confinement gaps between the continua of
1- & 2-string states.

Our TS3A approach offers a Bethe-string-based perspective, contributing to a
more integrated understanding of Heisenberg spin systems. This framework has
demonstrated its efficiency and validity in interpreting experimental
observations of the quasi-1D antiferromagnet $\rm YbAlO_{3}$ Yang  _et al._
(2023a). This intriguing consistency between theoretical predictions based on
TS3A and experimental results motivates its extended application to ladder and
two-dimensional systems featuring interchain coupling. This broadening of
scope not only enhances the versatility of the Bethe string picture but also
transcends its conventional one-dimensional limitations.

Acknowledgments.  We thank Yunfeng Jiang for helpful discussion. This work is
supported by National Natural Science Foundation of China No. 12274288 and the
Innovation Program for Quantum Science and Technology Grant No. 2021ZD0301900,
and the Natural Science Foundation of Shanghai with grant No. 20ZR1428400.

## References

  * Yang  _et al._ (2023a) J. Yang, T. Xie, S. E. Nikitin, J. Wu, and A. Podlesnyak, Phys. Rev. B 108, L020402 (2023a).
  * Jimbo and Miwa (1995) M. Jimbo and T. Miwa,  _Algebraic Analysis of Solvable Lattice Models_ (American Mathematical Society, 1995).
  * Franchini (2017) F. Franchini,  _An Introduction to Integrable Techniques for One-Dimensional Quantum Systems_, Vol. 940 (Springer, Cham, 2017).
  * Yang  _et al._ (2023b) J. Yang, X. Wang, and J. Wu, J. Phys. A: Math. Theor. 56, 013001 (2023b).
  * He  _et al._ (2017) F. He, Y. Jiang, Y.-C. Yu, H.-Q. Lin, and X.-W. Guan, Phys. Rev. B 96, 220401 (2017).
  * Gaudin (1967) M. Gaudin, Phys. Lett. A 24, 55 (1967).
  * Yang (1967) C. N. Yang, Phys. Rev. Lett. 19, 1312 (1967).
  * Guan  _et al._ (2013) X.-W. Guan, M. T. Batchelor, and C. Lee, Rev. Mod. Phys. 85, 1633 (2013).
  * Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963).
  * Lieb (1963) E. H. Lieb, Phys. Rev. 130, 1616 (1963).
  * Li  _et al._ (2023) R.-T. Li, S. Cheng, Y.-Y. Chen, and X.-W. Guan, J. Phys. A: Math. Theor. 56, 335204 (2023).
  * Guan and Andrei (2018) H. Guan and N. Andrei, arXiv.1803.04846 (2018).
  * Kitanine  _et al._ (1999) N. Kitanine, J. M. Maillet, and V. Terras, Nucl. Phys. B 554, 647 (1999).
  * Kitanine  _et al._ (2000) N. Kitanine, J. M. Maillet, and V. Terras, Nucl. Phys. B 567, 554 (2000).
  * Caux and Hagemans (2006) J.-S. Caux and R. Hagemans, J. Stat. Mech.: Theory Exp. 2006, P12013 (2006).
  * Caux  _et al._ (2008a) J.-S. Caux, J. Mossel, and I. P. Castillo, J. Stat. Mech.: Theory Exp. 2008, P08006 (2008a).
  * Caux  _et al._ (2008b) J.-S. Caux, J. Mossel, and I. P. Castillo, J. Stat. Mech.: Theory Exp. 2008, P08006 (2008b).
  * Castillo (2020) I. P. Castillo, arXiv:2005.10729 [cond-mat] (2020).
  * Takahashi (1971) M. Takahashi, Prog. Theor. Phys. 46, 401 (1971).
  * Kohno (2009) M. Kohno, Phys. Rev. Lett. 102, 037203 (2009).
  * Yang  _et al._ (2019) W. Yang, J. Wu, S. Xu, Z. Wang, and C. Wu, Phys. Rev. B 100, 184406 (2019).
  * Lake  _et al._ (2013) B. Lake, D. A. Tennant, J.-S. Caux, T. Barthel, U. Schollwöck, S. E. Nagler, and C. D. Frost, Phys. Rev. Lett. 111, 137205 (2013).
  * Wang  _et al._ (2018) Z. Wang, J. Wu, W. Yang, A. K. Bera, D. Kamenskyi, A. T. M. N. Islam, S. Xu, J. M. Law, B. Lake, C. Wu, and A. Loidl, Nature 554, 219 (2018).
  * Bera  _et al._ (2020) A. K. Bera, J. Wu, W. Yang, R. Bewley, M. Boehm, J. Xu, M. Bartkowiak, O. Prokhnenko, B. Klemke, A. T. M. N. Islam, J. M. Law, Z. Wang, and B. Lake, Nat. Phys. 16, 625 (2020).
  * Wang  _et al._ (2019) Z. Wang, M. Schmidt, A. Loidl, J. Wu, H. Zou, W. Yang, C. Dong, Y. Kohama, K. Kindo, D. I. Gorbunov, S. Niesen, O. Breunig, J. Engelmayer, and T. Lorenz, Phys. Rev. Lett. 123, 067202 (2019).
  * Steinigeweg  _et al._ (2016) R. Steinigeweg, J. Herbrych, X. Zotos, and W. Brenig, Phys. Rev. Lett. 116, 017202 (2016).
  * Batchelor  _et al._ (2005) M. T. Batchelor, X.-W. Guan, N. Oelkers, and Z. Tsuboi, Advances in Physics 56, 465 (2005).
  * Essler  _et al._ (1997) F. H. L. Essler, A. M. Tsvelik, and G. Delfino, Phys. Rev. B 56, 11001 (1997).
  * Starykh and Balents (2014) O. A. Starykh and L. Balents, Phys. Rev. B 89, 104407 (2014).
  * Bonča  _et al._ (1994) J. Bonča, J. P. Rodriguez, J. Ferrer, and K. S. Bedell, Phys. Rev. B 50, 3415 (1994).
  * Heidrich-Meisner  _et al._ (2003) F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W. Brenig, Phys. Rev. B 68, 134436 (2003).
  * Keselman  _et al._ (2020) A. Keselman, L. Balents, and O. A. Starykh, Phys. Rev. Lett. 125, 187201 (2020).
  * Zhou  _et al._ (2021) C. Zhou, Z. Yan, H.-Q. Wu, K. Sun, O. A. Starykh, and Z. Y. Meng, Phys. Rev. Lett. 126, 227201 (2021).
  * Takahashi and Suzuki (1972) M. Takahashi and M. Suzuki, Prog. Theor. Phys. 48, 2187 (1972).
  * Fujita  _et al._ (2003) T. Fujita, T. Kobayashi, and H. Takahashi, J. Phys. A: Math. Gen. 36, 1553 (2003).
  * Ilakovac  _et al._ (1999) A. Ilakovac, M. Kolanović, S. Pallua, and P. Prester, Phys. Rev. B 60, 7271 (1999).
  * Isler and Paranjape (1993) K. Isler and M. B. Paranjape, Phys. Lett. B 319, 209 (1993).
  * Essler  _et al._ (1992) F. H. L. Essler, V. E. Korepin, and K. Schoutens, J. Phys. A: Math. Gen. 25, 4115 (1992).
  * Hagemans and Caux (2007) R. Hagemans and J.-S. Caux, J. Phys. A: Math. Theor. 40, 14605 (2007).
  * Marboe and Volin (2017) C. Marboe and D. Volin, J. Phys. A: Math. Theor. 50, 204002 (2017).
  * Bajnok  _et al._ (2020) Z. Bajnok, E. Granet, J. L. Jacobsen, and R. I. Nepomechie, J. High Energ. Phys. 2020, 177 (2020).
  * Hou  _et al._ (2023) J. Hou, Y. Jiang, and R.-D. Zhu, arXiv.2303.07640 (2023).
  * Wu  _et al._ (2019) L. S. Wu, S. E. Nikitin, Z. Wang, W. Zhu, C. D. Batista, A. M. Tsvelik, A. M. Samarakoon, D. A. Tennant, M. Brando, L. Vasylechko, M. Frontzek, A. T. Savici, G. Sala, G. Ehlers, A. D. Christianson, M. D. Lumsden, and A. Podlesnyak, Nat. Commun. 10, 698 (2019).
  * Fan  _et al._ (2020) Y. Fan, J. Yang, W. Yu, J. Wu, and R. Yu, Phys. Rev. Res. 2, 013345 (2020).
  * Nikitin  _et al._ (2021) S. E. Nikitin, S. Nishimoto, Y. Fan, J. Wu, L. S. Wu, A. S. Sukhanov, M. Brando, N. S. Pavlovskii, J. Xu, L. Vasylechko, R. Yu, and A. Podlesnyak, Nat. Commun. 12, 3599 (2021).
  * Bera  _et al._ (2014) A. K. Bera, B. Lake, W.-D. Stein, and S. Zander, Phys. Rev. B 89, 094402 (2014).
  * Zou  _et al._ (2021) H. Zou, Y. Cui, X. Wang, Z. Zhang, J. Yang, G. Xu, A. Okutani, M. Hagiwara, M. Matsuda, G. Wang, G. Mussardo, K. Hódsági, M. Kormos, Z. He, S. Kimura, R. Yu, W. Yu, J. Ma, and J. Wu, Phys. Rev. Lett. 127, 077201 (2021).
  * Bethe (1931) H. Bethe, Z. Physik 71, 205 (1931).
  * Karabach  _et al._ (1997) M. Karabach, G. Müller, H. Gould, and J. Tobochnik, Comput. Phys. 11, 36 (1997).
  * Karbach  _et al._ (2000) M. Karbach, K. Hu, and G. Müller, arXiv:cond-mat/0008018 (2000).
  * Vladimirov (1986) A. A. Vladimirov, Theor. Math. Phys. 66, 102 (1986).
  * Takahashi (1999) M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models_ (Cambridge University Press, 1999).
  * (53) In the spin-1/2 XXX chain, we can check this invalidity by inserting the unphysical solutions into the wavefunction.
  * Deguchi and Giri (2015) T. Deguchi and P. R. Giri, J. Stat. Mech.: Theory Exp. 2015, P02004 (2015).
  * James  _et al._ (2018) A. J. A. James, R. M. Konik, P. Lecheminant, N. J. Robinson, and A. M. Tsvelik, Rep. Prog. Phys. 81, 046002 (2018).