USTC-ICTS/PCFT-23-43

Exact Spin Correlators of Integrable Quantum Circuits from Algebraic Geometry

Arthur Hutsalyuk1,2, Yunfeng Jiang3$\star$, Balázs Pozsgay2, Hefeng Xu4 and
Yang Zhang4,5

1 SISSA and INFN, Sezione di Trieste, via Bonomea 265, 34136 Trieste, Italy

2 MTA-ELTE “Momentum” Integrable Quantum Dynamics Research Group,

Department of Theoretical Physics, ELTE Eötvös Loránd University,

Pázmány Péter stny. 1A, Budapest 1117, Hungary

3 School of Physics and Shing-Tung Yau Center, Southeast University,

Nanjing 210096, Jiangsu, China

4 Interdisciplinary Center for Theoretical Study, University of Science and
Technology of China, Hefei, Anhui 230026, China

5 Peng Huanwu Center for Fundamental Theory, Hefei, Anhui 230026, China

⋆<EMAIL_ADDRESS>
## Abstract

We calculate the correlation functions of strings of spin operators for
integrable quantum circuits exactly. These observables can be used for
calibration of quantum simulation platforms. We use algebraic Bethe Ansatz, in
combination with computational algebraic geometry to obtain analytic results
for medium-size (around 10-20 qubits) quantum circuits. The results are
rational functions of the quantum circuit parameters. We obtain analytic
results for such correlation functions both in the real space and Fourier
space. In the real space, we analyze the short time and long time limit of the
correlation functions. In Fourier space, we obtain analytic results in
different parameter regimes, which exhibit qualitatively different behaviors.
Using these analytic results, one can easily generate numerical data to
arbitrary precision.

###### Contents

  1. 1 Introduction
  2. 2 The model and the observables
    1. 2.1 The model
    2. 2.2 The physical quantity
  3. 3 Bethe ansatz and algebraic geometry
    1. 3.1 Algebraic Bethe ansatz
    2. 3.2 Computational algebraic geometry
  4. 4 Exact results I. Real space
    1. 4.1 Results for small quantum numbers
    2. 4.2 Short time behavior
    3. 4.3 Long time behavior
    4. 4.4 A different initial state
  5. 5 Exact results II. Fourier space
    1. 5.1 Massive regime
    2. 5.2 Massless regime
  6. 6 Conclusions
  7. A Proof of overlap formula (70)
  8. B More on computational algebraic geometry
    1. B.1 Polynomial ring and ideal
    2. B.2 Gröbner basis and quotient ring
    3. B.3 Companion matrix
    4. B.4 Efficient implementations
  9. C Rational $Q$-system
  10. D Symmetrization
  11. E More analytic results in massive regime
  12. F More analytic results in massless regime

##  1 Introduction

Integrable models are special many body systems, which allow for the exact
computation of physical observables in various situations [1, 2]. These exact
results can then be compared to measurements in real world experiments. For
many decades, researches have focused on the physical quantities in
equilibrium, and in many cases spectacular agreement was found between theory
and experiment. Examples include the measurement of the exact free energy of
the 1D Bose gas [3], the computation of static correlation functions of the
Heisenberg spin chain [4], and the experimental observation of the
$E_{8}$-related spectrum of the perturbed Ising model [5, 6, 7, 8, 9].

However, in the last decade the focus of the research shifted towards non-
equilibrium situations. This was motivated by advances on the experimental
side: it became possible to realize well-known integrable models in various
types of experiments, enabling the study of non-equilibrium dynamics in these
systems. Examples for experiments with ultra-cold atoms include the famous
“Quantum Newton’s cradle” experiment [10], the observation of the Generalized
Gibbs Ensemble [11], and the confirmation of Generalized Hydrodynamics [12].

As an alternative to experiments with ultra-cold gases, it has also become
possible to simulate integrable models on quantum computers, as demonstrated
in the recent experiment of Google Quantum AI [13]. In a quantum computer the
time evolution is performed in discrete steps, via the succession of local
operations. It is possible to perform these steps in a homogeneous manner,
thus obtaining a so-called brickwork quantum circuit. Such circuits can be
integrable, and in such cases they are often related to well known statistical
physical models and quantum spin chains studied earlier in the literature [14,
15]

The integrable circuit studied in [13] is related to the six-vertex model and
the Heisenberg spin chain [16]. Even though these models are well known,
relatively few exact results have been published for the corresponding
brickwork circuits; for examples see [14, 17]. The Bethe Ansatz solution of
this specific circuit was considered in [18], where it was argued that a
certain correlation function of this circuit is a measurable quantity in
experiment. These correlation functions are crucial for error mitigation for
quantum simulation platforms [19]. However, the explicit computation of this
correlation function was not given in [18], except for eigenstates with only
two magnons.

In this work we revisit the problem of [18] and we compute closed form results
for the correlation function proposed there. We focus on spin chains with
intermediate sizes (10-20 qubits), which is appropriate for the quantum
circuits in real world experiments.

The correlation functions in question can be computed using the standard form
factor expansion, however such computation will involve summation over all
eigenstates of the finite volume system with given global quantum numbers.
Therefore, any theoretical method to compute these objects has to treat all
eigenstates, and in Bethe Anstaz this means to find all physical solutions to
the Bethe equations. This is a notoriously difficult task, see for example
[20]. Instead of solving the Bethe Ansatz equations numerically, we apply
recent developments in integrability, including the rational $Q$-system [21,
22, 23, 24, 25, 26, 27, 28, 29] and algebro-geometric approach [30, 31, 32,
33, 34], which allows us to calculate the final result _analytically_. We also
make use the algebraic Bethe ansatz, which allows us to use the full power of
integrability in the computation of correlation functions and obtain results
beyond two-magnons states.

The rest of the paper is structured as follows. In section 2, we introduce the
integrable quantum circuit model and the quantities that we would like to
compute. In section 3, we introduce the main tools for computations, namely
algebraic Bethe ansatz and computational algebraic geometry. The results are
presented in the following two sections. In section 4, we present results in
real space. We analyze the ‘short time’ limit and ‘long time’ limit where the
discrete ‘time’ here refers to the number of layers in the quantum circuit. We
found that, for short times the result is independent from the size of the
system if the system size is sufficiently large. In section 5 we present the
exact results in Fourier space 5. We computed the result explicitly for two
sets of parameters of the quantum circuit. They exhibit qualitatively
different behaviors. We conclude in section 6 and discuss future directions.
More technical details are delegated to the appendices.

##  2 The model and the observables

In this section, we introduce the model that we are going to study and the
physical quantities that we want to compute.

###  2.1 The model

We are interested in the dynamics of a one dimensional quantum cellular
automaton, where the time evolution is given by the successive application of
local two-site gates. The arrangement of these gates is depicted in figure 1;
such a system is often called a brickwork quantum circuit, and it can be
regarded as a special type of a Floquet system. For a formal definition of the
circuit, see below.

Figure 1: Graphical depiction of the quantum circuit in question. The vertical
direction is associated with the direction of time. Each rectangle stands for
the action of the unitary two-site gate $\check{R}$, given in the main text.

#### Evolution operator

The operator which evolves the whole system in time (depicted as the vertical
direction in figure 1) is defined as

$\displaystyle\mathcal{U}=\mathcal{U}_{\text{e}}\mathcal{U}_{\text{o}}$ (1)

where

$\displaystyle\mathcal{U}_{\text{o}}=$ $\displaystyle\,U_{12}U_{34}\ldots
U_{2L-1,2L}\,,$ (2) $\displaystyle\mathcal{U}_{\text{e}}=$
$\displaystyle\,U_{23}U_{45}\ldots U_{2L,1}\,,$

and $U_{j,j+1}$ stands for a two-site unitary operator (acting on sites $j$
and $j+1$) to be specified below.

#### Two-site Unitary

The operator $U$ that we consider was given in [18], and it reads

$\displaystyle U=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&\cos\alpha&\mathrm{i}\sin\alpha&0\\\ 0&\mathrm{i}\sin\alpha&\cos\alpha&0\\\
0&0&0&e^{2\mathrm{i}\phi}\end{array}\right)\,.$ (7)

Here $\alpha$ and $\phi$ are real numbers, such that $U$ is unitary. In the
spin chain language, this operator conserves the $z$ component of the spin,
which implies that the time evolution operator $\mathcal{U}$ will commute with
the global $S^{z}$ operator given by

$S_{z}=\frac{1}{2}(\sigma_{1}^{z}+\ldots+\sigma_{2L}^{z})\,.$ (8)

We will show that by a proper change of variables, this unitary gate is
equivalent to a certain evaluation of the so-called $\check{R}$-matrix known
in the integrability literature. To make this identification, we first
factorize out $\sigma^{z}$ operator and a phase factor

$\displaystyle
U=e^{\mathrm{i}\phi}e^{-\mathrm{i}\phi(\sigma_{1}^{z}+\sigma_{2}^{z})/2}\left(\begin{array}[]{cccc}1&0&0&0\\\
0&e^{-\mathrm{i}\phi}\cos\alpha&\mathrm{i}e^{-\mathrm{i}\phi}\sin\alpha&0\\\
0&\mathrm{i}e^{-\mathrm{i}\phi}\sin\alpha&e^{-\mathrm{i}\phi}\cos\alpha&0\\\
0&0&0&1\end{array}\right)\,.$ (13)

After multiplying all the $U$ operators, the factors in front simply yield a
global phase $e^{-\mathrm{i}\phi(2S^{z}+L)}$.

The so-called $\check{R}$-matrix of the XXZ spin chain is written in the
traditional notation as

$\displaystyle\check{R}(u)=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&c(u)&b(u)&0\\\ 0&b(u)&c(u)&0\\\ 0&0&0&1\end{array}\right)\,.$ (18)

The matrix elements of (18) read

$\displaystyle b(u)=\frac{\varphi(u)}{\varphi(u+\eta)},\qquad
c(u)=\frac{\varphi(\eta)}{\varphi(u+\eta)}\,,$ (19)

where $\varphi(u)=\sinh u$ and $\eta$ is related to the anisotropy parameter
$\Delta$ of the XXZ chain via $\Delta=\cosh\eta$. It is useful to consider
also the so-called rational limit corresponding to the $SU(2)$-symmetric XXX
chain, in which case $\varphi(u)=u$ and $\eta=\mathrm{i}$.

Now we can read off the identification between the two different
parametrizations given in (13) and (18). For the generic case we find

$\displaystyle
e^{-\mathrm{i}\phi}\cos\alpha=\frac{\sinh\eta}{\sinh(u+\eta)},\qquad\mathrm{i}e^{-\mathrm{i}\phi}\sin\alpha=\frac{\sinh
u}{\sinh(u+\eta)}\,.$ (20)

We can easily solve for $\alpha$ and $\phi$ in terms of $u$ and $\eta$ as
follows

$\displaystyle\tan\alpha=-\mathrm{i}\frac{\sinh u}{\sinh\eta},\qquad
e^{-2\mathrm{i}\phi}=-\frac{\sinh(u-\eta)}{\sinh(u+\eta)}\,.$ (21)

The rational limit is obtained only when $\alpha=\phi$.

In [18], the parameters $\alpha$ and $\phi$ are real numbers, but for our
intermediate computations we can work with generic complex $u$ and $\eta$. In
this parametrization, $u$ has the interpretation of the spectral parameter.

From now on, we shall work with the $u$, $\eta$ parametrization and write the
two-site unitary as $\check{R}_{ij}(u)$ as in the integrability literature
where $i$ and $j$ denote the sites that it acts on.

The most important property of $\check{R}_{ij}$ is that it satisfies quantum
Yang-Baxter equation

$\displaystyle\check{R}_{ij}(u)\check{R}_{ik}(u+v)\check{R}_{jk}(v)=\check{R}_{jk}(v)\check{R}_{ik}(u+v)\check{R}_{ij}(u)\,,$
(22)

which underlies the integrability of the model.

#### Remarks

This model is well known in the literature of integrable models, and it made
its appearance in several different contexts, which we briefly summarize
below.

  1. 1.

It first appeared in classical integrable lattice models by putting the
6-vertex model on the diagonal lattice [16]. In this context, $\check{R}_{ij}$
encodes the local Boltzmann weights and $\mathcal{U}$ is known as the
diagonal-to-diagonal transfer matrix, which can be used to compute lattice
partition functions.

  2. 2.

The above model also serve as a lattice discretization of the quantum sine-
Gordon model, which was use by Destri and de Vega to compute the finite size
energy spectrum [35, 36, 37]. In this context, $\mathcal{U}$ is related to the
evolution in the light cone directions.

  3. 3.

The model appeared in [14] as the integrable Trotterization of the XXZ spin
chain. and further studied in [15].

We want to emphasis that although the underlying model is equivalent, the
contexts are quite different. The interesting physical quantities depend on
the context. We will now turn to the quantities that we are interested in the
current context.

###  2.2 The physical quantity

In statistical mechanics and condensed matter physics the research is often
focused on the properties of the ground states or the first few low-lying
excited states. These computations are typically done in the thermodynamic
limit. However, in the current context we are considering a quantum circuit
with finite qubits and discrete time evolution. Therefore, there is no notion
of a ground state. Furthermore, we will not consider the thermodynamic limit,
because we intend to compute exact results with a fixed finite number of
qubits, or spins111In what follows we will often use the language of spin
chains. Therefore we shall use qubits and spins interchangeably..

The quantity that can be measured conveniently for quantum computers are
correlation functions of local spin operators acting on the pseudovacuum state
[18]. Let us denote the basis states of each qubit by $|0\rangle$ and
$|1\rangle$. The pseudovacuum state is defined by

$\displaystyle|\Omega\rangle=|0\rangle_{1}\otimes|0\rangle_{2}\otimes\ldots\otimes|0\rangle_{2L}\equiv|00\ldots
0\rangle\,.$ (23)

We consider the following string of local spin operators

$\displaystyle\hat{\Sigma}_{k}^{(N)}=\sigma_{k+1}^{x}\sigma_{k+2}^{x}\ldots\sigma_{k+N}^{x}$
(24)

which acts on $N$ sites starting from site-$k$.

Acting on the pseudo-vacuum, this operator creates a domain wall state. For
concreteness, let us define

$\displaystyle|\text{DW}_{N}\rangle\equiv\hat{\Sigma}_{0}^{(N)}|\Omega\rangle=|\underbrace{11\ldots
1}_{N}0\ldots 0\rangle\,.$ (25)

We are interested in the dynamical correlation function

${\left\langle\Omega\right|}\hat{\Sigma}_{k}^{(N)}(0)\hat{\Sigma}_{0}^{(N)}(n){\left|\Omega\right\rangle},$
(26)

where time evolution of the operators is understood in the Heisenberg picture:

$\hat{\Sigma}_{0}^{(N)}(n)\equiv U^{n}\hat{\Sigma}_{0}^{(N)}(0)U^{-n}$ (27)

Apart from trivial phases, the pseudo-vaccum is a fixed point of the time
evolution operator, therefore the correlation can also be expressed as

$\displaystyle\mathcal{D}^{(N)}(n,k)=\langle\Omega|\hat{\Sigma}_{k}^{(N)}\mathcal{U}^{n}\hat{\Sigma}_{0}^{(N)}|\Omega\rangle\,.$
(28)

For $k=0$, this quantity is expressed further as

${\left\langle\text{DW}_{N}\right|}\mathcal{U}^{n}{\left|\text{DW}_{N}\right\rangle},$
(29)

which is the discrete version of the so-called Loschmidt echo or return
amplitude (fidelity) in quantum quench dynamics.

It is also useful to consider the Fourier transform of $\mathcal{D}^{(N)}$
defined by222Notice that in the summand we have $\mathcal{D}^{(N)}(n,2k)$
instead of $\mathcal{D}^{(N)}(n,k)$. This definition is more natural because
the system is invariant under a 2-site translation.

$\displaystyle\widehat{\mathcal{D}}^{(N)}(\omega,Q)=\sum_{n=0}^{\infty}e^{-n(\mathrm{i}\omega+\gamma)}\sum_{k=1}^{L}e^{\mathrm{i}Qk}\mathcal{D}^{(N)}(n,2k)$
(30)

where the real parameter $\gamma\ll 1$ is introduced to provide proper
convergence for the time Fourier transform. Here $Q$ is called the quasi-
momentum and it satisfies the quantization condition $e^{\mathrm{i}QL}=1$.

Our goal is computing the correlation functions (28) and (30) analytically for
given quantum numbers.

##  3 Bethe ansatz and algebraic geometry

Analytic computation of quantities like (28) and (30) are usually quite hard.
Thanks to the integrability of the model, we can obtain analytic results for
these quantities using Bethe ansatz and computational algebraic geometry. In
this section, we discuss how to compute observables of the integrable quantum
circuits using these powerful theoretical tools.

###  3.1 Algebraic Bethe ansatz

#### Transfer matrices

In order to compute time evolution $\mathcal{U}$ analytically, we diagonalize
it by Bethe ansatz and work with its eigenvalues. In the integrability
literature, $\mathcal{U}$ is identified with the so-called diagonal-to-
diagonal transfer matrix. Its explicit diagonalization is achieved by a well-
known relation between the diagonal-to-diagonal transfer matrix and the row-
to-row transfer matrix, which we will briefly recall now. We start with the
$R$-matrix of the XXZ spin chain

$\displaystyle R_{ij}(u)=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&b(u)&c(u)&0\\\ 0&c(u)&b(u)&0\\\ 0&0&0&1\end{array}\right)_{ij}$ (35)

This $R$-matrix is related to the $\check{R}$-matrix by
$\check{R}_{ij}(u)=P_{ij}R_{ij}(u)$ where $P_{ij}$ is the permutation operator

$\displaystyle P_{ij}=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&0&1&0\\\
0&1&0&0\\\ 0&0&0&1\end{array}\right)_{ij}\,.$ (40)

Following the standard procedure of the algebraic Bethe ansatz (see _e.g_
[38]), we can define the following monodromy matrices

$\displaystyle M_{0}(u|\\{\theta_{j}\\})=$
$\displaystyle\,R_{01}(u-\theta_{1})R_{02}(u-\theta_{2})\ldots
R_{0,2L}(u-\theta_{2L})\,,$ (41)
$\displaystyle\widehat{M}_{0}(u|\\{\theta_{j}\\})=$
$\displaystyle\,R_{0,2L}(u+\theta_{2L})\ldots
R_{0,2}(u+\theta_{2})R_{0,1}(u+\theta_{1})\,.$

where ‘$0$’ denotes an auxiliary space and we have introduced a set of
parameters $\\{\theta_{j}\\}$ called the inhomogeneities. The row-to-row
transfer matrix is defined as

$\displaystyle\mathcal{T}(u|\\{\theta_{j}\\})=\text{tr}_{0}M_{0}(u|\\{\theta_{j}\\})\,,\qquad\widehat{\mathcal{T}}(u|\\{\theta_{j}\\})=\text{tr}_{0}\widehat{M}_{0}(u|\\{\theta_{j}\\})\,.$
(42)

Our evolution operator $\mathcal{U}$ defined in (1), (2) is related to the
row-to-row transfer matrix (42) by [39]

$\displaystyle\mathcal{U}(u)=\mathcal{T}\left(\tfrac{u}{2}|\theta_{j}(\tfrac{u}{2})\right)\widehat{\mathcal{T}}\left(\tfrac{u}{2}|\theta_{j}(\tfrac{u}{2})\right)$
(43)

where $\theta_{j}(\tfrac{u}{2})$ stands for the following choice of the
inhomogeneities

$\displaystyle\theta_{j}(\tfrac{u}{2})=\left\\{-\tfrac{u}{2},\tfrac{u}{2},\ldots,-\tfrac{u}{2},\tfrac{u}{2}\right\\}\,.$
(44)

The relation (43) can be proven by using the so-called regularity condition
$R_{ij}(0)=P_{ij}$, which is valid for our $R$-matrix (35).

The transfer matrix $\widehat{\mathcal{T}}$ is related to the transfer matrix
$\mathcal{T}$ by

$\displaystyle\mathcal{T}(-u-\eta|\\{\theta_{j}\\})=\widehat{\mathcal{T}}(u|\\{\theta_{j}\\})\,,$
(45)

which can be proven straightforwardly by using the crossing relation of the
$R$-matrix

$\displaystyle
R_{12}(-u-\eta)=-\sigma_{1}^{y}R_{12}^{t_{1}}(u)\sigma_{1}^{y}\,.$ (46)

The row-to-row transfer matrix $\mathcal{T}$ can be diagonalized by Bethe
ansatz. We denote each eigenstate of $\mathcal{T}$ by $|\mathbf{u}_{K}\rangle$
where $\mathbf{u}_{K}=\\{u_{1},u_{2},\ldots,u_{K}\\}$ are the Bethe roots that
are rapidities of magnons. We have

$\displaystyle\mathcal{T}(u|\\{\theta_{j}\\})|\mathbf{u}_{K}\rangle=t(u|\\{\theta_{j}\\})|\mathbf{u}_{K}\rangle$
(47)

where the eigenvalue is given by

$\displaystyle
t(u|\\{\theta_{j}\\})=\prod_{k=1}^{K}\frac{\varphi(u-u_{k}-\eta)}{\varphi(u-u_{k})}+\prod_{j=1}^{2L}\frac{\varphi(u-\theta_{j})}{\varphi(u-\theta_{j}+\eta)}\prod_{k=1}^{K}\frac{\varphi(u-u_{k}+\eta)}{\varphi(u-u_{k})}\,.$
(48)

The Bethe roots satisfy Bethe ansatz equations (BAE)

$\displaystyle\prod_{a=1}^{2L}\frac{\varphi(u_{k}-\theta_{a}+\eta)}{\varphi(u_{k}-\theta_{a})}=\prod_{j\neq
k}^{K}\frac{\varphi(u_{k}-u_{j}+\eta)}{\varphi(u_{k}-u_{j}-\eta)},\qquad
k=1,\ldots,K\,.$ (49)

Using the relation (43) and (45), we obtain the eigenvalue of $\mathcal{U}$

$\displaystyle\mathcal{U}|\mathbf{u}_{K}\rangle=\tau(\mathbf{u}_{K})|\mathbf{u}_{K}\rangle$
(50)

where

$\displaystyle\tau(\mathbf{u}_{K})=\prod_{k=1}^{K}\frac{\varphi(u_{k}-\tfrac{u}{2}+\eta)}{\varphi(u_{k}-\tfrac{u}{2})}\frac{\varphi(u_{k}+\tfrac{u}{2})}{\varphi(u_{k}+\tfrac{u}{2}+\eta)}\,.$
(51)

For the choice of inhomogeneities (44), the BAE (49) becomes

$\displaystyle\left(\frac{\varphi(u_{k}-\tfrac{u}{2}+\eta)}{\varphi(u_{k}-\tfrac{u}{2})}\right)^{L}\left(\frac{\varphi(u_{k}+\tfrac{u}{2}+\eta)}{\varphi(u_{k}+\tfrac{u}{2})}\right)^{L}=\prod_{j\neq
k}^{K}\frac{\varphi(u_{k}-u_{j}+\eta)}{\varphi(u_{k}-u_{j}-\eta)}\,,\qquad
k=1,\ldots,K.$ (52)

#### Shift operators

It is also useful to define the shift operator

$\displaystyle V=P_{12}P_{23}\ldots P_{2L-1,2L}$ (53)

which shifts all qubits towards right by one site. We have

$\displaystyle\hat{\Sigma}^{(N)}_{k}|\Omega\rangle=V^{k}\hat{\Sigma}_{0}^{(N)}|\Omega\rangle=V^{k}|\text{DW}_{N}\rangle\,.$
(54)

Due to our choice of inhomogeneities (44) which is alternating, the one-site
shift operator $V$ is not compatible with most Bethe states. In our context,
it is more natural to define the two-site shift operator $\mathcal{V}=V^{2}$.
The Bethe states are eigenstates of $\mathcal{V}$:

$\displaystyle\mathcal{V}|\mathbf{u}_{K}\rangle=\mu(\mathbf{u}_{K})|\mathbf{u}_{K}\rangle$
(55)

where the eigenvalue $\mu(\mathbf{u}_{K})$ is given by

$\displaystyle\mu(\mathbf{u}_{K})=\prod_{k=1}^{K}\frac{\varphi(u_{k}-\frac{u}{2})}{\varphi(u_{k}-\tfrac{u}{2}+\eta)}\frac{\varphi(u_{k}+\frac{u}{2})}{\varphi(u_{k}+\tfrac{u}{2}+\eta)}$
(56)

Due to the periodic boundary condition, $\mu(\mathbf{u}_{K})$ satisfies the
following quantization condition

$\displaystyle\mu(\mathbf{u}_{K})^{L}=1\,.$ (57)

#### Overlaps

Inserting a resolution of identity in (28), we obtain

$\displaystyle\mathcal{D}^{(N)}(n,2k)=\sum_{\text{sol}_{N}}\frac{\tau(\mathbf{u}_{N})^{n}}{\mu(\mathbf{u}_{N})^{k}}\frac{\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle\langle\mathbf{u}_{K}|\text{DW}_{N}\rangle}{\langle\mathbf{u}_{N}|\mathbf{u}_{N}\rangle}$
(58)

where ‘solN’ means that we take the sum over all physical solutions of the
Bethe ansatz equations (52) with $N$ Bethe roots. The reason that we only need
to sum over the $N$-magnon sector is that $|\text{DW}_{N}\rangle$ only has
non-zero overlap with $N$-magnon states, due to $S^{z}$ conservation. To write
the result in terms of the Bethe roots, we also need analytic expressions for
the overlaps $\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle$,
$\langle\mathbf{u}_{N}|\text{DW}_{N}\rangle$ and the norm
$\langle\mathbf{u}_{N}|\mathbf{u}_{N}\rangle$. Fortunately, these quantities
have been computed [40, 41, 42, 43]. We summarize the results below.

The overlap is given by

  1. 1.

Even $N$

$\displaystyle\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle=$
$\displaystyle\,\prod_{k=1}^{N}\varphi(u_{k}-\tfrac{u}{2}+\eta)^{N/2}\varphi(u_{k}+\tfrac{u}{2}+\eta)^{N/2}$
(59) $\displaystyle\times\frac{\det
H^{\text{even}}}{\varphi(u)^{N^{2}/4}\prod_{j>k}\varphi(u_{j}-u_{k})}$

where $H^{\text{even}}$ is an $N\times N$ matrix whose matrix elements are
given by

$\displaystyle
H_{ab}^{\text{even}}=\left\\{\begin{array}[]{ll}g_{b}(u_{a},\tfrac{u}{2}),&1\leq
b\leq\tfrac{N}{2}\\\ g_{b-\frac{N}{2}}(u_{a},-\tfrac{u}{2}),&\frac{N}{2}<b\leq
N\end{array}\right.$ (62)

with

$\displaystyle
g_{n}(u,\theta)=\frac{\varphi^{\prime}(u-\theta)}{\varphi(u-\theta)}-\frac{\varphi^{\prime}(u-\theta+\eta)}{\varphi(u-\theta+\eta)}\,.$
(63)

  2. 2.

Odd $N$

$\displaystyle\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle=$
$\displaystyle\,\prod_{k=1}^{N}\varphi(u_{k}-\tfrac{u}{2}+\eta)^{\frac{N+1}{2}}\varphi(u_{k}+\tfrac{u}{2}+\eta)^{\frac{N-1}{2}}$
(64) $\displaystyle\times\frac{\det
H^{\text{odd}}}{\varphi(u)^{(N^{2}-1)/4}\prod_{j>k}^{N}\varphi(u_{j}-u_{k})}$

where $H^{\text{odd}}$ is an $N\times N$ matrix whose matrix elements are
given by

$\displaystyle
H_{ab}^{\text{odd}}=\left\\{\begin{array}[]{ll}g_{b}(u_{a},\tfrac{u}{2}),&1\leq
b\leq\tfrac{N+1}{2}\\\
g_{b-\frac{N+1}{2}}(u_{a},-\tfrac{u}{2}),&\frac{N+1}{2}<b\leq
N\end{array}\right.$ (67)

and the function $g_{a}(u,\theta)$ is given in (63).

The norm of the on-shell Bethe state is given by

$\displaystyle\langle\mathbf{u}_{N}|\mathbf{u}_{N}\rangle=\varphi(\eta)^{N}\prod_{k=1}^{N}\left(\frac{\varphi(u_{k}-\tfrac{u}{2})}{\varphi(u_{k}-\tfrac{u}{2}+\eta)}\right)^{2L}\left(\frac{\varphi(u_{k}+\tfrac{u}{2})}{\varphi(u_{k}+\tfrac{u}{2}+\eta)}\right)^{2L}\prod_{j\neq
k}^{N}\frac{\varphi(u_{j}-u_{k}+\eta)}{\varphi(u_{j}-u_{k})}\times\det G$ (68)

where $G$ is the Gaudin matrix whose elements are

$\displaystyle G_{ab}=-\frac{\partial}{\partial
u_{b}}\ln\left[\left(\frac{\varphi(u_{a}-\tfrac{u}{2}+\eta)}{\varphi(u_{a}-\tfrac{u}{2})}\right)^{L}\left(\frac{\varphi(u_{a}+\tfrac{u}{2}+\eta)}{\varphi(u_{a}+\tfrac{u}{2})}\right)^{L}\right]-\frac{\partial}{\partial
u_{b}}\ln\left[\prod_{k=1}^{N}\frac{\varphi(u_{a}-u_{k}-\eta)}{\varphi(u_{a}-u_{k}+\eta)}\right]\,.$
(69)

Finally we also need the overlap $\langle\mathbf{u}_{N}|\text{DW}_{N}\rangle$.
This is the complex conjugate of $\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle$.
However, it is not convenient to work with complex conjugate in the
computational algebraic geometry method. For even $N$, we have the following
relation

$\displaystyle\langle\mathbf{u}_{N}|\text{DW}_{N}\rangle=\prod_{k=1}^{N}\left(\frac{\varphi(u_{k}-\tfrac{u}{2})}{\varphi(u_{k}-\tfrac{u}{2}+\eta)}\right)^{N/2}\left(\frac{\varphi(u_{k}+\tfrac{u}{2})}{\varphi(u_{k}+\tfrac{u}{2}+\eta)}\right)^{N/2}\,\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle\,.$
(70)

The proof is given in appendix A. Defining the following quantity

$\displaystyle
w(\mathbf{u}_{N})=\frac{\langle\text{DW}_{n}|\mathbf{u}_{N}\rangle\langle\mathbf{u}_{N}|\text{DW}_{N}\rangle}{\langle\mathbf{u}_{N}|\mathbf{u}_{N}\rangle}$
(71)

where the analytic expressions for the three overlaps have been given in (59),
(70) and (68). Therefore our task is to compute the following quantity

$\displaystyle{\mathcal{D}^{(N)}(n,2k)=\sum_{\text{sol}_{N}}\frac{\tau(\mathbf{u}_{N})^{n}}{\mu(\mathbf{u}_{N})^{k}}w(\mathbf{u}_{N})}\,.$
(72)

for arbitrary $n$. The Fourier transform of this quantity can be computed
readily

$\displaystyle\widehat{\mathcal{D}}^{(N)}(\omega,Q)=$
$\displaystyle\,\sum_{k=1}^{L}e^{\mathrm{i}kQ}\sum_{n=0}^{\infty}e^{-n(\mathrm{i}\omega+\gamma)}{\mathcal{D}}^{(N)}(n,2k)$
(73) $\displaystyle=$
$\displaystyle\,\sum_{\text{sol}_{N}}\left(\sum_{k=1}^{L}\frac{e^{\mathrm{i}kQ}}{\mu(\mathbf{u}_{N})^{k}}\right)\frac{w(\mathbf{u}_{N})}{1-e^{-\mathrm{i}\omega-\gamma}\tau(\mathbf{u}_{N})}.$

Since $\mu(\mathbf{u}_{N})$ are quantized, we can write it as

$\displaystyle\mu(\mathbf{u}_{N})=e^{\mathrm{i}Q_{0}},\qquad Q_{0}=\frac{\pi
m}{L},\qquad m=0,1,\ldots,L-1\,.$ (74)

If we fix the total quasi-momentum $Q=Q_{0}$, (73) simplifies to

$\displaystyle{\widehat{\mathcal{D}}^{(N)}(\omega,Q_{0})=L\,\sum_{\text{sol}_{N,Q_{0}}}\frac{w(\mathbf{u}_{N})}{1-e^{-\mathrm{i}\omega-\gamma}\tau(\mathbf{u}_{N})}}$
(75)

where now the summation is over the physical solutions of the Bethe ansatz
equations with total quasi-momentum $Q_{0}$. This can be achieved by adding
the additional constraint

$\displaystyle
e^{\mathrm{i}Q_{0}}=\prod_{k=1}^{N}\frac{\varphi(u_{k}-\tfrac{u}{2})}{\varphi(u_{k}-\tfrac{u}{2}+\eta)}\frac{\varphi(u_{k}+\tfrac{u}{2})}{\varphi(u_{k}+\tfrac{u}{2}+\eta)}\,$
(76)

to the BAE. This will in general reduce the number of allowed solutions.

###  3.2 Computational algebraic geometry

In order to compute (72) and (75), one of the main difficulties is summing
over solutions of Bethe ansatz equations. In principle, this can only be done
numerically by solving BAE and plugging in the solutions into (72) and (75).
However, there are drawbacks with this approach. The BAE and the physical
quantities depend on the parameters $\eta$ and $u$. In order to solve the BAE
numerically, we need to fix these parameters to certain values. On the other
hand, computations of small lattice size shows that the final results are
rational functions in these parameters333More precisely, for the isotropic
case $\eta=0$, the result are rational functions in $u$. For $\eta\neq 0$, the
results are rational functions of $e^{\eta}$ and $e^{u}$.. In many situations,
it is much more desirable to have analytic expressions so that one can perform
different analysis. Therefore, we would like to obtain results (72) and (75)
that are as analytical as possible, with explicit $u$ or $\eta$ (or both)
dependence. Numerical solution cannot achieve this goal.

To this end, we apply the algebro-geometric method which allows us to compute
the sum over solutions of BAE in (72) and (75) analytically. This approach is
completely analytic and exact, which avoids solving the equations explicitly.
The efficiency of the computation depends on its implementations. With proper
implementation, we can reach obtain analytic results for $L\sim 10$
(corresponding to 20 qubits444Note that the number of qubits in the brickwork
model is $2L$.) with various magnon numbers. Our approach nicely covers the
range of medium quantum numbers where brute force computations become
infeasible while thermodynamic approximation is not always applicable.

#### Rational $Q$-system

To perform the sum in (72) and (75) analytically, there is one more
difficulty. It is well-known that among solutions of BAE there exist ones that
do not correspond to eigenstates of the transfer matrix, such solutions are
called non-physical solutions. They have to be discarded while performing the
sum. Deciding which solutions are physical is a non-trivial task and is
related to the completeness problem of Bethe ansatz. In order to sum over only
the physical solutions, we apply the method of rational $Q$-system. The
rational $Q$-system is a non-trivial reformulation of the BAE, which is more
efficient than solving BAE555We compare efficiency for obtaining _all_
physical solutions using the two different formulations. In terms of finding
the ground state energy, the original BAE has certain advantages.. More
importantly, it leads to only physical solutions. This was first observed in
[21] and was proven more rigorously in [22, 23, 28]. For XXX and XXZ spin
chains, the rational $Q$-systems have been constructed in [21] and [23]
respectively.

We give a brief introduction to the rational $Q$-system in appendix C. Here it
is sufficient to mention that, in rational $Q$-system, we are solving a set of
algebraic equations for the coefficient of Baxter’s $Q$-polynomials. For a
state with $N$ Bethe roots, Baxter’s $Q$-polynomials for XXX spin chain
(corresponding to $\eta=0$) as

$\displaystyle
Q_{\text{XXX}}(u)=\prod_{j=1}^{N}(u-u_{j})=u^{N}+\sum_{k=0}^{N-1}c_{k}\,u^{k}\,.$
(77)

The rational $Q$-system will give a set of algebraic equations for
$\\{c_{k}\\}$. For XXZ chain (corresponding to $\eta\neq 0$), Baxter’s
polynomial is defined by

$\displaystyle\widetilde{Q}_{\text{XXZ}}(u)=\prod_{j=1}^{N}\sinh(u-u_{j})$
(78)

in terms of additive variables $u$ and $u_{k}$. For the XXZ case, it is more
convenient to use the multiplicative variables

$\displaystyle x=e^{2u},\qquad x_{k}=e^{2u_{k}},\qquad
q=e^{\eta},\qquad\xi=e^{2\theta}\,.$ (79)

The $Q$-function in (78) becomes

$\displaystyle\widetilde{Q}_{\text{XXZ}}(u)\propto
Q_{\text{XXZ}}(x)=\prod_{j=1}^{N}(x-x_{j})=x^{N}+\sum_{k=0}^{N-1}c_{k}\,x^{k}\,.$
(80)

In the framework of rational $Q$-system, one first derives a set of equations
for the coefficients $\\{c_{k}\\}$ in (77) or (80). These algebraic equations
are called the zero remainder conditions. One then either solve these
equations explicitly, or manipulate these equations analytically using
computational algebraic geometric approaches.

#### Symmetrization

Once the coefficients $\\{c_{k}\\}$ are known, we have all the
$Q$-polynomials. To find the Bethe roots, we simply need to find the zeros of
$Q$-polynomials. However, since our aim is computing physical quantities
analytically, we can also avoid this procedure. All physical quantities are
symmetric functions of Bethe roots. More precisely, the quantities we are
considering are _symmetric rational functions_ of Bethe roots. At the same
time, from (77) it is clear that $\\{c_{k}\\}$ are essentially the elementary
symmetric polynomials of $\\{u_{k}\\}$, _i.e._

$\displaystyle c_{N-1}=$ $\displaystyle\,-(u_{1}+u_{2}+\ldots+u_{N})\,,$ (81)
$\displaystyle c_{N-2}=$
$\displaystyle\,u_{1}u_{2}+u_{1}u_{3}+\ldots+u_{N-1}u_{N}\,,$
$\displaystyle\cdots$ $\displaystyle c_{0}=$
$\displaystyle\,(-1)^{N}\,u_{1}u_{2}\ldots u_{N}\,.$

Therefore it is clear that all the physical quantities can also be written in
terms of $\\{c_{k}\\}$.

In the algebraic geometrical computation, we will work with expressions in
terms of $\\{c_{k}\\}$. Since the overlaps are originally given in terms of
the Bethe roots $\\{u_{k}\\}$, it is necessary to perform a change of
variables from $\\{u_{k}\\}$ to $\\{c_{k}\\}$. When the expressions are
relatively simple with a few variables, it is straightforward to perform the
calculation666For example, using the function SymmetricReduction of
Mathematica.. However, for more involved expressions, _e.g._ explicit analytic
formula of the Gaudin determinant, such a procedure can be quite tedious. At
the moment we do not have an analytic formula for Gaudin-like determinants
directly written in terms of $\\{c_{k}\\}$. Nevertheless, we can perform the
change of variables with the help of some tricks using the Groebner basis,
which accelerates the computation significantly. We will discuss the
symmetrization procedure in more details in Appendix D.

#### Computational algebraic geometry

Now we discuss the main procedure for the algebro-geometric computations. More
details can be found in appendix B.

Suppose we want to compute the following sum

$\displaystyle
S=\sum_{\text{sol}}\frac{P(x_{1},\ldots,x_{n})}{Q(x_{1},\ldots,x_{n})}$ (82)

where $P(x_{1},\ldots,x_{n})$ and $Q(x_{1},\ldots,x_{n})$ are polynomials in
$x_{1},\ldots,x_{n}$. Here $\\{x_{k}\\}$ are solutions of a set of polynomial
equations

$\displaystyle f_{k}(x_{1},\ldots,x_{n})=0,\qquad k=1,\ldots,N.$ (83)

We assume that the system (83) has finitely many solutions, which is the case
for Bethe equations for $\\{u_{k}\\}$, or equivalently the zero remainder
conditions for $\\{c_{k}\\}$. To perform the sum in (82), usually we first
solve the equations (83) and find all possible solutions, and then plug the
solutions to the summand in (82) and finally perform the sum. In the AG
computation, the strategy is different. In order to explain the main idea, let
us introduce a few notions in algebraic geometry (or commutative algebra). All
polynomials in $n$ variables $\\{x_{1},\ldots,x_{n}\\}$ with complex
coefficients form a polynomial ring, which we denote by
$\mathbb{C}[x_{1},\ldots,x_{n}]$. The polynomials in (83) which define the
polynomial equations generate an ideal

$\displaystyle\mathcal{I}=\langle f_{1},f_{2},\ldots,f_{N}\rangle\,.$ (84)

In other words, the ideal $\mathcal{I}$ consists of all polynomials of the
form $\sum_{i}a_{i}f_{i}$ with $a_{i}\in\mathbb{C}[x_{1},\ldots,x_{n}]$. One
important fact of ideal $\mathcal{I}$ is that the same ideal can be generated
by different set of polynomials. These polynomials are called bases of the
ideal. The most important kind of bases for us is the so-called Gröbner basis
(see details in Appendix B.2), which we will denote by $g_{j}$ and we have

$\displaystyle\mathcal{I}=\langle f_{1},f_{2},\ldots,f_{N}\rangle=\langle
g_{1},g_{2},\ldots,g_{N^{\prime}}\rangle$ (85)

where in general $N^{\prime}\neq N$. The Gröbner basis has the nice feature
that polynomial reductions with respect to it has well defined and unique
remainders.

The quotient space

$\displaystyle V_{f}=\mathbb{C}[x_{1},\ldots,x_{n}]/\mathcal{I}$ (86)

is called the quotient ring. The quotient ring is a finite dimensional vector
space and the dimension $\dim V_{f}=N_{s}$ equals to the number of solutions
of the set of equations (83). Since $V_{f}$ is an $N_{s}$-dimensional vector
space, we can choose a basis which we denote by $e_{1},\ldots,e_{N_{s}}$.
There is a canonical choice of such a basis, which can be constructed from the
Gröbner basis of $\mathcal{I}$. After introducing the canonical basis $e_{j}$,
any polynomial $P\in\mathbb{C}[x_{1},\ldots,x_{n}]$ can be mapped to a finite
dimensional matrix $\mathbf{P}$ in $V_{f}$ by

$\displaystyle P\,e_{i}\sim\sum_{j=1}^{N_{s}}[\mathbf{P}]_{ij}\,e_{j}$ (87)

where ‘$\sim$’ means performing polynomial reduction with respect to the
Gröbner basis of the ideal $\mathcal{I}$ and finding the remainder. The
remainder can be written as a linear combination of the basis $e_{j}$ with
coefficients that are matrix elements $[\mathbf{P}]_{ij}$. The $N_{s}\times
N_{s}$ dimensional matrix $\mathbf{P}$ is called _companion matrix_ of the
polynomial $P$. It can be seen as the ‘representation’ of $P$ in $V_{f}$.

To summarize, given the system of equations (83) and a polynomial
$P(x_{1},\ldots,x_{n})$, we construct its companion matrix $\mathbf{P}$ by the
following steps

  1. 1.

Compute the Gröbner basis of the equations (83);

  2. 2.

Using the Gröbner basis, construct the canonical basis of the quotient ring;

  3. 3.

Using both the Gröbner basis and quotient ring basis, construct the companion
matrix.

All the above computations can be automatized. For the computation of Gröbner
basis, we use the package Singular. The quotient ring and companion matrices
can be computed using our own Mathematica package BeF2, which is currently
private. In certain computation for larger quantum numbers, similar procedures
are also implemented with Singular and GPI-Space framework to boost efficiency
of certain computation of companion matrix both over real and Fourier space.

Now coming back to our sum (82). For the two polynomials
$P(x_{1},\ldots,x_{n})$ and $Q(x_{1},\ldots,x_{n})$ in the summand, we can
construct their companion matrices

$\displaystyle P(x_{1},\ldots,x_{n})\mapsto\mathbf{P},\qquad
Q(x_{1},\ldots,x_{n})\mapsto\mathbf{Q}$ (88)

The crucial results that we need from computational algebraic geometry is
that, the sum (82) can be computed analytically by

$\displaystyle S=\text{tr}\left[\mathbf{P}\cdot\mathbf{Q}^{-1}\right]$ (89)

where $\mathbf{Q}^{-1}$ is the inverse matrix of $\mathbf{Q}$.

#### Dealing with free parameters

Before ending the section, let us make a remark. One of the main challenges in
the current context for the AG computation is that our equations have free
parameters $u$ and $\eta$. In principle, the algorithms also work for such
cases. In practice, however, the procedure usually involve bulky intermediate
expressions which slows down the computation. Therefore we need a better way
to implement the computation. The strategy we take is to perform the
computation for different rational values of the parameters and then fit the
final result.

##  4 Exact results I. Real space

In this section, we present analytic results for the correlation function

$\displaystyle\mathcal{D}^{(N)}(n,k)=\langle\Omega|\hat{\Sigma}_{k}^{(N)}\mathcal{U}^{n}\hat{\Sigma}_{0}^{(N)}|\Omega\rangle\,.$
(90)

We present the results for $k=0$, in which case $\mathcal{D}^{(N)}(n,0)$ is
essentially a discrete version of the Loschmidt echo. Results for $k\neq 0$
can be obtained easily because the shift in $k$ only leads to a global phase
factor $e^{\mathrm{i}kQ_{0}}$. Written in terms of rapidities, it is given by

$\displaystyle\mathcal{D}^{(N)}(n,0)=\sum_{\text{sol}_{N}}{\tau(\mathbf{u}_{N})^{n}}\,w(\mathbf{u}_{N})\,.$
(91)

Before presenting the result, let us discuss the general structure. From the
definition (90), it is clear that the final result is given by certain matrix
elements of $\mathcal{U}^{n}$. From the definition of $\mathcal{U}$ (43),
matrix elements are polynomials of $b(u)$ and $c(u)$ where $b(u)$ and $c(u)$
are given in (19). The result depends on three global quantum numbers: the
size of the circuit, which we denote by $2L$; the length of the string of
operators $\hat{\Sigma}^{(N)}$ which is denoted by $N$; and the evolution time
denoted by $n$. Therefore we expect the results take the following form

$\displaystyle\mathcal{D}_{L}^{(N)}(n)=\sum_{l,k}\alpha_{l,k}\left(\frac{\varphi(u)}{\varphi(u+\eta)}\right)^{l}\left(\frac{\varphi(\eta)}{\varphi(u+\eta)}\right)^{k},\qquad
l,k\geq 0$ (92)

where $\alpha_{l,k}$ are non-negative integers. In the rational limit,
$\varphi(u)\to u$ and $\eta\to\mathrm{i}$, we find

$\displaystyle\mathcal{D}^{(N)}_{L}(n)=\sum_{l,k}\alpha_{l,k}\left(\frac{u}{u+\mathrm{i}}\right)^{l}\left(\frac{\mathrm{i}}{u+\mathrm{i}}\right)^{k},\qquad
l,k\geq 0$ (93)

Therefore, in the isotropic limit, the result is a rational function of $u$.
The AG computation leads to such rational functions. Below we will discuss the
results in the isotropic case (93). Once we obtain the results in the rational
limit using AG computation, we can solve for the coefficients $\alpha_{l,k}$.
Plugging $\alpha_{l,k}$ back to (92) gives the result for generic $\eta$.

###  4.1 Results for small quantum numbers

For small values of $L,N,n$, we can compute $\mathcal{D}_{L}^{(N)}(n)$ by
brute force, without using Bethe Ansatz. These results can serve as benchmarks
for the Bethe Ansatz and AG computation. They already exhibit some interesting
feature. We list the results for $2L=4,6,8$ below for small $N$ and $n$.

#### Results for $L=2$

For $N=1$, we have

$\displaystyle\mathcal{D}^{(1)}_{2}(1)=$ $\displaystyle\,c^{2}\,,$ (94)
$\displaystyle\mathcal{D}^{(1)}_{2}(2)=$ $\displaystyle\,(b^{2}+c^{2})^{2}\,,$
$\displaystyle\mathcal{D}^{(1)}_{2}(3)=$
$\displaystyle\,(3b^{2}c+c^{3})^{2}\,,$
$\displaystyle\mathcal{D}^{(1)}_{2}(4)=$
$\displaystyle\,(b^{4}+6b^{2}c^{2}+c^{4})^{2}\,,$
$\displaystyle\mathcal{D}^{(1)}_{2}(5)=$
$\displaystyle\,(5b^{4}c+10b^{2}c^{3}+c^{5})^{2}\,.$

where $b,c$ are short for $b(u)$ and $c(u)$. For $N=2$, we obtain

$\displaystyle\mathcal{D}^{(2)}_{2}(1)=$ $\displaystyle\,c^{2}\,,$ (95)
$\displaystyle\mathcal{D}^{(2)}_{2}(2)=$
$\displaystyle\,2b^{4}c^{2}+b^{4}+2b^{2}c^{4}+c^{4}\,,$
$\displaystyle\mathcal{D}^{(2)}_{2}(3)=$
$\displaystyle\,2b^{8}c^{2}+6b^{6}c^{4}+4b^{6}c^{2}+6b^{4}c^{6}+16b^{4}c^{4}+3b^{4}c^{2}+2b^{2}c^{8}+4b^{2}c^{6}+c^{6}\,,$
$\displaystyle\mathcal{D}^{(2)}_{2}(4)=$
$\displaystyle\,b^{8}+6b^{8}c^{2}+4b^{10}c^{2}+2b^{12}c^{2}+6b^{4}c^{4}+42b^{6}c^{4}+40b^{8}c^{4}+10b^{10}c^{4}+42b^{4}c^{6}$
$\displaystyle\,+72b^{6}c^{6}+20b^{8}c^{6}+c^{8}+6b^{2}c^{8}+40b^{4}c^{8}+20b^{6}c^{8}+4b^{2}c^{10}+10b^{4}c^{10}+2b^{2}c^{12}$

#### Results for $L=3$

For $N=1$, we have

$\displaystyle\mathcal{D}^{(1)}_{3}(1)=$ $\displaystyle\,c^{2}\,,$ (96)
$\displaystyle\mathcal{D}^{(1)}_{3}(2)=$ $\displaystyle\,2b^{2}c^{2}+c^{4}\,,$
$\displaystyle\mathcal{D}^{(1)}_{3}(3)=$
$\displaystyle\,b^{6}+3b^{4}c^{2}+6b^{2}c^{4}+c^{6}\,,$
$\displaystyle\mathcal{D}^{(1)}_{3}(4)=$
$\displaystyle\,12b^{6}c^{2}+18b^{4}c^{4}+12b^{2}c^{6}+c^{8}\,.$

For $N=2$, we have

$\displaystyle\mathcal{D}^{(2)}_{3}(1)=$ $\displaystyle\,c^{2}\,,$ (97)
$\displaystyle\mathcal{D}^{(2)}_{3}(2)=$
$\displaystyle\,b^{4}c^{2}+c^{4}+2b^{2}c^{4}\,,$
$\displaystyle\mathcal{D}^{(2)}_{3}(3)=$
$\displaystyle\,b^{8}+4b^{8}c^{2}+4b^{4}c^{4}+12b^{6}c^{4}+c^{6}+4b^{2}c^{6}+9b^{4}c^{6}+2b^{2}c^{8}\,$

#### Results for $L=4$

For $N=1$, we have

$\displaystyle\mathcal{D}^{(1)}_{4}(1)=$ $\displaystyle\,c^{2}\,,$ (98)
$\displaystyle\mathcal{D}^{(1)}_{4}(2)=$ $\displaystyle\,2b^{2}c^{2}+c^{4}\,,$
$\displaystyle\mathcal{D}^{(1)}_{4}(3)=$
$\displaystyle\,3b^{4}c^{2}+6b^{2}c^{4}+c^{6}$

For $N=2$, we have

$\displaystyle\mathcal{D}^{(2)}_{4}(1)=$ $\displaystyle\,c^{2}\,,$ (99)
$\displaystyle\mathcal{D}^{(2)}_{4}(2)=$
$\displaystyle\,b^{4}c^{2}+c^{4}+2b^{2}c^{4}\,,$

###  4.2 Short time behavior

At short times, the system does not evolve much and cannot ‘feel’ the size the
of system. Therefore we expect for small $n$, the result should be independent
of the length $L$. This is indeed what we observed from our analytic results.
Let us illustrate this using $N=2$ examples. Larger magnon numbers exhibit
similar behaviors. For $N=2$, we find that

$\displaystyle\mathcal{D}_{2}^{(2)}(1)=\mathcal{D}_{3}^{(2)}(1)=\ldots=\mathcal{D}_{L}^{(2)}(1)=-\frac{1}{(u+\mathrm{i})^{2}},\qquad
L\geq 2.$ (100)

Namely the result is the same for $n=1$ for all $L\geq 2$. Similarly, we find
for $n=2$

$\displaystyle\mathcal{D}_{3}^{(2)}(2)=\mathcal{D}_{4}^{(2)}(2)=\ldots\mathcal{D}_{L}^{(2)}(2)=\frac{-u^{4}+3u^{2}+2\mathrm{i}\,u-1}{(u+\mathrm{i})^{6}},\qquad
L\geq 3$ (101)

and

$\displaystyle\mathcal{D}_{4}^{(2)}(3)=$
$\displaystyle\,\mathcal{D}_{5}^{(2)}(3)=\ldots=\mathcal{D}_{L}^{(2)}(3),\quad
L\geq 4$ (102) $\displaystyle=$
$\displaystyle\,\frac{-u^{8}+12u^{6}+8\mathrm{i}\,u^{5}-18u^{4}-12\mathrm{i}\,u^{3}+12u^{2}+4\mathrm{i}\,u-1}{(u+\mathrm{i})^{10}}\,.$

In general, the pattern is that for any $n\geq 1$, the quantity
$\mathcal{D}^{(2)}_{L}(n)$ is independent of $L$ for $L\geq n$. In other
words, this indicates that in $n$ discrete steps the system evolve at most to
size $2n$. If the system is large enough, this spread does not touch the
boundary of the system and therefore the result is independent of $L$.

###  4.3 Long time behavior

As one can envisage, the quantity $\mathcal{D}^{(N)}_{L}(n)$ will be more
involved as the discrete time $n$ grows. For large $n$, the analytic result we
obtain is a rather complicated rational function of $u$. The brute force
computation becomes infeasible, but the result can be computed by the AG
method. In this Subsection and in later parts of this work we present results
obatined by the AG methods.

In order to gain some idea what the result looks like, we present an example
with $L=8$, $N=2$ and $n=10$

$\displaystyle\mathcal{D}^{(2)}_{8}(10)=\frac{P(u)}{(u+\mathrm{i})^{38}}$
(103)

where $P(u)$ is the following polynomial

$\displaystyle P(u)=$
$\displaystyle\,-1+18\mathrm{i}\,u+243u^{2}-1776\mathrm{i}\,u^{3}-12345u^{4}+57576\mathrm{i}\,u^{5}+276288u^{6}$
(104)
$\displaystyle\,-939576\mathrm{i}\,u^{7}-3394500u^{8}+8772432\mathrm{i}\,u^{9}+25101788u^{10}$
$\displaystyle\,-50778808\mathrm{i}\,u^{11}-119366908u^{12}+196227896\mathrm{i}\,u^{13}+382479184u^{14}$
$\displaystyle\,-523609032\mathrm{i}\,u^{15}-837663006u^{16}+961297284\mathrm{i}\,u^{17}+1247817562u^{18}$
$\displaystyle\,-1194362808\mathrm{i}\,u^{19}-1250264310u^{20}+986339608\mathrm{i}\,u^{21}+831213408u^{22}$
$\displaystyle\,-531759240\mathrm{i}\,u^{23}-361426868u^{24}+183920576\mathrm{i}\,u^{25}+101357548u^{26}$
$\displaystyle\,-40023464\mathrm{i}\,u^{27}-18065212u^{28}+5290568\mathrm{i}\,u^{29}+1977488u^{30}$
$\displaystyle\,-391416\mathrm{i}\,u^{31}-121065u^{32}+13530\mathrm{i}\,u^{33}+3435u^{34}-120\mathrm{i}\,u^{35}-25u^{36}\,.$

Although $n=10$ is not a very large number, the polynomial $P(u)$ is already
rather bulky. For larger values of $n$, say around $n=100$, the result
involves huge polynomials in $u$. Presenting the results explicitly is not
very illuminating. In order to gain more intuitions about the long time
behavior of the Loschmidt echo, we instead plot the zeros of the corresponding
$\mathcal{D}^{(N)}_{L}(n)$ as a function of $u$. These zeros are akin to the
Lee-Yang zeros of thermal partition functions. It is known that the Lee-Yang
zeros encode interesting physics of the model such as the phase structure,
here we also want to investigate how the zeros of $\mathcal{D}^{(N)}_{L}(n)$
changes with $n$.

In figure 2(a) and 2(b), we plot the zeros of $\mathcal{D}^{(2)}_{4}$ for
$n=6$ and $n=10$. We can see that the zeros both distribute symmetrically
about the imaginary axis. The two distribution of zeros look somewhat similar,
but does not exhibit clear pattern.

(a) The zeros for $L=4$, $M=2$, $n=6$

(b) The zeros for $L=4$, $M=2$, $n=10$

However, as $n$ increases, the zeros will start to condense and a clear
pattern emerges, as is shown in figure 3.

Figure 3: The zeros of $\mathcal{D}_{4}^{(2)}(500)$. There are 1998 zeros in
total and most of them condense on three curves as can be seen clearly, one
along the real axis and the rest two in the lower half plane.

This phenomena, namely the condensation of zeros of the
$\mathcal{D}_{L}^{(M)}(n)$ for large $n$ is also manifested for other quantum
numbers, with different condensation curves. In figure 4, we show the zeros of
$\mathcal{D}_{4}^{(4)}(500)$. Comparing the condensation curves of
$\mathcal{D}_{4}^{(2)}(500)$ and $\mathcal{D}_{4}^{(4)}(500)$, we find that
the condensation curves have more complicated structure for larger $M$.

Figure 4: The zeros of $\mathcal{D}_{4}^{(4)}(500)$. There are 3992 zeros in
total and they clearly condense on a few condensation curves, one along the
real axis and two others in the lower half plane.

###  4.4 A different initial state

The distribution of zeros that we have seen in the previous subsection is very
similar to the zeros of the partition function of the lightcone 6-vertex
model. In fact, we can identify the two quantities as follows. Let us consider
the following state

$\displaystyle\widetilde{\Sigma}|\Omega\rangle=\prod_{j=1}^{L}\left(|01\rangle-|10\rangle\right)^{\otimes}=\prod_{j=1}^{L}\left(\sigma_{2j-1}^{x}-\sigma_{2j}^{x}\right)|00\ldots
0\rangle\,,$ (105)

which is the so-called dimer state. We can define a similar correlation
function to (90)

$\displaystyle\widetilde{\mathcal{D}}_{L}(n)=\langle\Omega|\widetilde{\Sigma}\,\mathcal{U}^{n}\,\widetilde{\Sigma}|\Omega\rangle\,.$
(106)

Interestingly, the above correlation function $\widetilde{\mathcal{D}}_{L}(n)$
coincides with the cylinder partition function of the 6-vertex model (up to
simple normalization factor) in statistical mechanics. For the rational case,
it has been studied in [32] and the zeros of the partition functions have been
studied in the partial thermodynamic limit. We have the following
identification $\widetilde{\mathcal{D}}_{L}(n)\propto Z_{L,n}$ where
$Z_{L,n}(u)$ is the partition function with size $L$ and $n$ in two
directions. $Z_{L,n}(u)$ is a polynomial in $u$. The partial thermodynamic
limit correspond to the limit with fixed finite $L$ and large $n$, which
corresponds precisely to the long time behavior of the correlation function
(106). One finds that the zeros of $Z_{L,n}(u)$ or equivalently
$\widetilde{\mathcal{D}}_{L}(n)$ also condense on certain curves for large
$n$. We compute the zeros of $\widetilde{\mathcal{D}}_{L}(n)$ for $L=4$,
$n=500$ in figure 5.

Figure 5: The zeros of $\widetilde{\mathcal{D}}_{8}(500)$, there are 3996
zeros in total, the zeros aslo condensates on a few condensation curves. The
condensation curve is similar to the one of $\mathcal{D}_{4}^{(4)}(500)$. An
obvious difference is that it also has condensation curves on the upper half
plane.

In statistical mechanics, the condensation curves separate different phases of
the model. One can even plot the condensation curves using numerical
approaches and find a nice agreement with the distribution of the zeros [32].

Interestingly, we find that the condensation curves of
$\widetilde{\mathcal{D}}_{4}(500)$ and $\mathcal{D}_{4}^{(4)}(500)$ are very
similar in the lower half plane. This can be seen clearly by putting the two
figures on top of each other, as is shown in figure 6 .

Figure 6: Comparison of the zeros of $\widetilde{\mathcal{D}}_{4}(500)$ and
$\mathcal{D}_{4}^{(4)}(500)$. We see that the condensation curves in the lower
half plane almost coincide with each other.

However, there are no condensation curves on the upper half plane for
$\mathcal{D}_{4}^{(4)}(500)$777Note that this does not mean that there are no
roots with positive imaginary parts. It means that the majority of such roots
are distributed very near the real axis. For example, there are only 50 roots
of $\mathcal{D}_{4}^{(4)}(500)$ with imaginary part larger $0.15$, which
constitute around $1\%$ of the roots. As a comparison, there are 1704 roots
with imaginary parts which are smaller than $-0.15$, which constitute $42\%$
of the roots.. The similarity between the condensation curves on the lower
half plane can be understood as follows. The appearance of condensation curve
in the long time limit is a consequence of the Beraha-Kahane-Weiss theorem
[44], which applies o the expressions of the form

$\displaystyle\mathcal{Z}_{n}(u)=\sum_{i}\alpha(u)\Lambda_{i}(u)^{n}$ (107)

where we shall refer to $\Lambda_{i}(u)$ as eigenvalues and $\alpha_{i}(u)$ as
the corresponding multiplicities. Comparing to (91), we see that in the
current context, $\Lambda_{i}(u)=\tau(\mathbf{u}_{N})$ are the eigenvalues of
the evolution operator $\mathcal{U}$ and the corresponding multiplicity
$\alpha_{i}(u)=w(\mathbf{u}_{N})$. We can arrange the eigenvalues according to
the norm as follows

$\displaystyle|\Lambda_{1}(u)|\geq|\Lambda_{2}(u)|\geq\ldots$ (108)

We say an eigenvalue is dominant if its norm is maximal. The BKW theorem
states that in the $n\to\infty$ limit, the zeros of $\mathcal{Z}_{n}(u)$ will
either form isolated points or curves. An isolated accumulation point occurs
at $u=z_{0}$ if there is a unique dominant eigenvalue, namely
$|\Lambda_{1}(z_{0})|>|\Lambda_{2}(z_{0})|$. A curve of accumulation points
occur when there are at least two dominant eigenvalues
$|\Lambda_{2}(u)|=|\Lambda_{2}(u)|$ and the relative phase
$\phi(u)\in\mathbb{R}$ defined by
$\Lambda_{2}(u)=e^{\mathrm{i}\phi(u)}\Lambda_{1}(u)$ varies along the curve.
Therefore the condensation curve will appear where there are at least two
eigenvalues dominant. Both $\widetilde{\mathcal{D}}_{4}(n)$ and
$\mathcal{D}_{4}^{(4)}(n)$ are computed in the $L=8,M=4$ sector and share the
same set of eigenvalues $\Lambda_{i}(u)$. Therefore a partial overlap of the
condensation curves is expected from the BKW theorem. On the other hand, why
there are no condensation curves for $\mathcal{D}_{4}^{(4)}(n)$ on the upper
half plane is unclear for us so far. It would be an interesting question to
have a better understanding on this point.

##  5 Exact results II. Fourier space

In this section, we consider the result of exact correlation functions in
Fourier space

$\displaystyle{\widehat{\mathcal{D}}^{(N)}(\omega,Q_{0})=L\,\sum_{\text{sol}_{N,Q_{0}}}\frac{w(\mathbf{u}_{N})}{1-e^{-\mathrm{i}\omega-\gamma}\tau(\mathbf{u}_{N})}}.$
(109)

The sum is now performed over the $N$ magnon sector with definite quasi-
momentum $e^{\mathrm{i}Q_{0}}$. For concreteness, we take the value $Q_{0}=0$,
which is the zero momentum sector for the given magnon number $N$. Other
momentum sectors can be computed similarly. The AG computation procedure is
similar, but the dimension of the quotient ring is smaller because we focus on
the zero momentum sector. Let us denote the companion matrices of
$w(\mathbf{u}_{N})$ and $\tau(\mathbf{u}_{N})$ by

$\displaystyle
w(\mathbf{u}_{N})\mapsto\mathbf{M}_{w},\qquad\tau(\mathbf{u}_{N})\mapsto\mathbf{M}_{\tau}\,.$
(110)

The analytic result is given by

$\displaystyle\widehat{\mathcal{D}}^{(N)}_{L}(\omega,0)=L\,\text{tr}\,\left[\mathbf{M}_{w}\cdot(1-a\,\mathbf{M}_{\tau})^{-1}\right],\qquad
a=e^{-\mathrm{i}\omega-\gamma}.$ (111)

The quantity $\widehat{\mathcal{D}}^{(N)}_{L}(\omega,0)$ is the one studied in
[18] numerically. We consider the generic XXZ $\check{R}$-matrix. In general,
the result depends on $\xi=e^{u}$, $q=e^{\eta}$ and
$a=e^{-\mathrm{i}\omega-\gamma}$ and is a rational function of these three
variables. Even for relatively small quantum numbers, the explicit analytic
result is rather involved. As an example, the analytic result for $L=4$, $M=2$
reads

$\displaystyle\widehat{\mathcal{D}}^{(2)}_{2}(\omega,0)=\frac{N_{2,2}(\xi,q,a)}{D_{2,2}(\xi,q,a)}$
(112)

where

$\displaystyle N_{2,2}(\xi,q,a)$
$\displaystyle=\,-a^{2}(q-\xi)^{2}(\xi+q)^{2}\left(\xi^{2}+q^{2}\left(\xi^{4}+\xi^{2}\left(q^{2}-4\right)+1\right)\right)$
(113)
$\displaystyle\,+a\left(\xi^{2}+q^{2}\left(\xi^{4}+\xi^{2}\left(q^{2}-4\right)+1\right)\right)\left(\xi^{2}+q^{2}\left(\xi^{4}\left(q^{2}+1\right)+\xi^{2}\left(q^{2}-6\right)+1\right)+1\right)$
$\displaystyle\,-\left(\xi^{2}q^{2}-1\right)^{4}$

and

$\displaystyle D_{2,2}(\xi,q,a)$
$\displaystyle=\,a^{3}\left(q^{2}-\xi^{2}\right)^{4}$ (114)
$\displaystyle\,-a^{2}\left(2\xi^{4}+\xi^{2}+\left(\xi^{2}+2\right)q^{4}+\left(\xi^{4}-8\xi^{2}+1\right)q^{2}\right)\left(\xi^{2}+q^{2}\left(\xi^{4}+\xi^{2}\left(q^{2}-4\right)+1\right)\right)$
$\displaystyle\,+a\left(\xi^{2}+q^{2}\left(\xi^{4}+\xi^{2}\left(q^{2}-4\right)+1\right)\right)\left(\xi^{2}+q^{2}\left(\xi^{4}\left(2q^{2}+1\right)+\xi^{2}\left(q^{2}-8\right)+1\right)+2\right)$
$\displaystyle\,-\left(\xi^{2}q^{2}-1\right)^{4}.$

In what follows, we will fix $\xi$ and $q$ to certain values and present the
result as a rational function of $a$. We distinguish two different regimes,
corresponding to $\Delta<1$ and $\Delta>1$ where
$\Delta=\tfrac{1}{2}(q+q^{-1})$. Borrowing the terminology from the XXZ spin
chain, we call $\Delta<1$ and $\Delta>1$ as the massless and massive regimes,
respectively. We will see that the results in the two different regimes are
qualitatively different.

###  5.1 Massive regime

We first consider the massive regime which is simpler. We take $\xi=9$ and
$q=2$. According to (21), this corresponds to the regime where $\alpha$ and
$\phi$ are complex numbers. we find that

$\displaystyle\widehat{\mathcal{D}}^{(2)}_{2}(\omega,0)=-\frac{156104641a^{2}-3440094482a+10884540241}{35153041a^{3}-1005425523a^{2}+6186972723a-10884540241}\,.$
(115)

We can plot the absolute value $|\widehat{\mathcal{D}}^{(2)}_{L}(\omega,0)|$
as a function of $\omega$, which is given in figure 7.

Figure 7: Plot of $|\widehat{\mathcal{D}}^{(2)}_{2}(\omega,0)|$ with
parameters $\xi=9$ and $q=2$. We plot in the range $-2\pi\leq\omega\leq 2\pi$.

Taking the same parameters $\xi=9$ and $q=2$, we present a few more results

$\displaystyle\widehat{\mathcal{D}}^{(2)}_{L}(\omega,0)=\frac{\hat{P}_{L}^{(2)}(a)}{\hat{Q}_{L}^{(2)}(a)}.$
(116)

For $L=3$, $M=2$, we have

$\displaystyle\hat{P}_{3}^{(2)}(a)=$
$\displaystyle\,25626566889a^{3}+2808199326933a^{2}$
$\displaystyle\,+88100949869467a-1135573198803289\,,$
$\displaystyle\hat{Q}_{3}^{(2)}(a)=$
$\displaystyle\,-208422380089a^{4}+4607080618756a^{3}+5564104859466a^{2}$
$\displaystyle\,+96035779705156a-1135573198803289.$

For $L=4$, $M=2$, we have

$\displaystyle\hat{P}_{4}^{(2)}(a)=$
$\displaystyle\,607759660339524a^{4}+418511232917558704a^{3}-7464677848782143656a^{2}$
$\displaystyle-38414127109876466896a+473892865031793352324$
$\displaystyle\hat{Q}_{4}^{(2)}(a)=$
$\displaystyle\,-4942945166190724a^{5}+114584571922332820a^{4}$
$\displaystyle\,+936172980724473560a^{3}-14807141945943110360a^{2}$
$\displaystyle-41725458557586857620a+473892865031793352324.$

For $L=5$, $M=2$, we have

$\displaystyle\hat{P}_{5}^{(2)}(a)=$
$\displaystyle\,900851756538259449a^{5}+29597987586042784755a^{4}$
$\displaystyle\,+2984695966367841261690a^{3}+13961140703024682671910a^{2}$
$\displaystyle\,+1044915399498613698674845a-12360192178975492163652649$
$\displaystyle\hat{Q}_{5}^{(2)}(a)=$
$\displaystyle\,-7326680472586200649a^{6}+177733069582320607894$
$\displaystyle a^{5}+329949279946254002265a^{4}$
$\displaystyle\,+5941087402154920804980a^{3}+22556459386698417302265a^{2}$
$\displaystyle\,+1131282374150658037135894a-12360192178975492163652649\,.$

Taking $a=e^{\mathrm{i}\omega}$, the corresponding plot is given in figure 8.

Figure 8: Plot of $|\widehat{\mathcal{D}}^{(2)}_{L}(\omega,0)|$ for $L=3\sim
10$. The curves converge to a steady shape as $L$ increases.

For higher quantum numbers, we can obtain similar exact results. The results
for $M=2$ and $L=3\sim 5$ can be found in the appendix E. We also obtained
analytic results for 4-magnon states for $L=4,5,6,7,8$. These results except
for $L=4,M=4$ are too bulky to be presented in the paper. Interested readers
can find this data in the attached file of arXiv submission. We present the
plot for the Fourier space correlation functions.

For 2-magnon states, as $L$ increases, the analytic results become more
complicated as expected. The number of solutions that we need to take into
account increases. This is given in table 1.

$L$ | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---
$\mathcal{N}_{s}$ | 5 | 8 | 9 | 12 | 13 | 16 | 17 | 20
Table 1: Number of solutions that contribute for $L=3\sim 10$, $M=2$.

However, as we can see from the plot, the shape of the curve quickly
converges. For large enough $L$, the difference between
$|\widehat{\mathcal{D}}^{(2)}_{L}(\omega,0)|$ and
$|\widehat{\mathcal{D}}^{(2)}_{L+1}(\omega,0)|$ become smaller. Their
difference is plotted in figure 9.

Figure 9: Plot of
$|\widehat{\mathcal{D}}^{(2)}_{L+1}(\omega,0)|-|\widehat{\mathcal{D}}^{(2)}_{L}(\omega,0)|$
for $L=6\sim 9$.

As we can see, the difference
$|\widehat{\mathcal{D}}^{(2)}_{10}(\omega,0)|-|\widehat{\mathcal{D}}^{(2)}_{9}(\omega,0)|$
is of order $10^{-6}$. This implies that for relatively large $L$, the result
is already a good approximation for the infinite length $L\to\infty$.

For 4-magnon states, we have the similar observation. The computation is
significantly more involved. The number of Bethe roots which contributes are
given in Table 2.

$L$ | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---
$\mathcal{N}_{s}$ | 20 | 42 | 85 | 143 | 232 | 340
Table 2: Number of Bethe roots that contributes for $L=4\sim 9$ and $M=4$.

The curves converge faster than the 2-magnon states. The plot for $L=4\sim 8$
is given in figure 11.

Figure 10: Plot of $|\widehat{\mathcal{D}}^{(4)}_{L}(\omega,0)|$ for $L=4\sim
8$. The curves converge so quickly and we cannot distinguish the curves for
$L\geq 5$ in the figure.

The differences are

$\displaystyle|\widehat{\mathcal{D}}^{(2)}_{7}(\omega,0)|-|\widehat{\mathcal{D}}^{(2)}_{6}(\omega,0)|\sim
10^{-5},\qquad|\widehat{\mathcal{D}}^{(2)}_{8}(\omega,0)|-|\widehat{\mathcal{D}}^{(2)}_{7}(\omega,0)|\sim
10^{-8}.$ (117)

They are plotted in figure 11.

Figure 11: Plot of
$|\widehat{\mathcal{D}}^{(4)}_{L+1}(\omega,0)|-|\widehat{\mathcal{D}}^{(4)}_{L}(\omega,0)|$
for $L=5\sim 7$.

###  5.2 Massless regime

In the massless regime, we take $\xi=9$,
$q=\tfrac{3}{5}+\tfrac{4\mathrm{i}}{5}$. According to (21), this corresponds
to the regime where $\alpha$ and $\phi$ are real numbers. The behaviors of
$|\widehat{\mathcal{D}}_{L}^{(M)}(\omega,0)|$ are rather different from the
massive regime. For fixed $M$, with increasing $L$, the curves do not
obviously converge. For $M=2$, we plot the
$|\widehat{\mathcal{D}}_{L}^{(2)}(\omega,0)|$ with $L=4,5,6,7,8,9$ in figure
12. Their comparison is given in figure 13

Figure 12: The plot of $|\widehat{\mathcal{D}}_{L}^{(2)}(\omega,0)|$ for
$L=4,5,6,7,8,9$. We take the range $-\pi\leq\omega<\pi$ with $\gamma=0.01$.
Figure 13: Comparison of $|\widehat{\mathcal{D}}_{L}^{(2)}(\omega,0)|$ with
$L=4,5,6,7,8,9$. the curves do not obviously converge as $L$ increases.

The analytic expressions for $\widehat{\mathcal{D}}_{L}^{(2)}(\omega,0)$ can
be found in the appendix. Similarly, we have obtained the analytic results for
$M=4$ with $L=4,5,6,7,8,9$. Their comparison is given in figure 15. The
analytic expression of the simplest one, namely
$\widehat{\mathcal{D}}_{4}^{(4)}(\omega,0)$ are two long to be presented
explicitly on the paper. Interested readers can download them from the
auxiliary file. We plot the results of
$\widehat{\mathcal{D}}_{L}^{(4)}(\omega,0)$ in figure 14 .

Figure 14: The plot of $|\widehat{\mathcal{D}}_{L}^{(4)}(\omega,0)|$ for
$L=4,5,6,7,8$. We take the range $-\pi\leq\omega<\pi$ and $\gamma=0.01$.
Figure 15: Comparison of $|\widehat{\mathcal{D}}_{L}^{(4)}(\omega,0)|$ with
$L=4,5,6,7,8,9$. The curves do not obviously converge as $L$ increases.

##  6 Conclusions

In this paper, we developed an analytic approach to compute the correlation
functions for a string of operators in an integrable quantum circuit. The
analytic result could be used to test real world implementations of quantum
circuits. Our method is based on Bethe ansatz solution of the model and
computational algebraic geometry, which allows us to sum over all solutions of
BAE analytically. The results we obtain are exact and analytic, without any
approximations. We obtained results for medium size systems with 10-20 qubits,
which is precisely the typical size of quantum circuits at the present stage.

We obtain results in both real space and Fourier space. In the real space, the
results of the correlation functions are rational functions of parameters
$\xi$ and $q$ (in the isotropic limit it is a a function of $u$). Using the
analytic results, we investigate both the short time and long time behaviors
of the correlation functions.

In the short time limit, we find that the result does not depend on the size
of the system if the size is larger than $2n$ where $n$ is the discrete time.
This is in agreement with the light cone propagation of information in the
brickwork quantum circuits.

In the long time limit for fixed system size, the result is equivalent to
partial thermodynamic limit of the cylinder partition function of the 6-vertex
model. Viewing the final result as a function of the spectral parameter, we
investigate its distribution of zeros. We found that the zeros condense on a
few condensation curves, similar to the zeros of the partition function of the
6-vertex model. This is a consequence of the BKW theorem. We also compared the
results to another choice of the string of operators, which corresponds to the
so-called dimer state in the XXZ spin chain. For $L=8$, $M=4$, we found that
the condensation curves for the two results overlap in the lower half plane,
but is quite different on the upper half plane. This interesting phenomena
deserves further investigation.

In Fourier space, the results are rational functions of $\xi$, $q$ and
$a=e^{\mathrm{i}\omega}$. We can obtain analytic results for generic fixed
$\xi$ and $q$, which is a rational function in $a$. The plot of such result
should be comparable to experimental data. We obtain these results for the
length $L=4\sim 9$ (corresponding to $8-18$ qubits) with fixed parameters
$\xi$ and $q$. We obtained results for both two- and four-magnon states. We
investigated two sets of parameters $(\xi,q)=(9,2),(2,3/5+4\mathrm{i}/5)$.
Using the terminology from the XXZ spin chain, they are in the massive and
massless regime, respectively. Using the parametrization of the quantum
circuits, they correspond to the case where $\alpha,\phi$ being complex and
real. We observed that the results are qualitatively different. For the
massive regime, the results are analytic for real $\omega$. For fixed magnon
number $M$, the numerical results for quickly converge with increasing $L$.
Therefore $L=9$ is already a very good approximation for the result in the
thermodynamic limit. In the massless regime, the result as a function of
$\omega$ have poles, therefore we need to set small non-zero $\gamma$ in
$a=e^{\mathrm{i}\omega+\gamma}$ to regularize the result. We find that for
fixed magnon numbers, the results do not converge quickly and the results are
rather different for different $L$. Therefore it is questionable whether in
such a case we can apply the thermodynamic approximation for finite $L$
result.

There are several interesting questions that one can investigate following the
current work. In the current paper, we mainly computed the correlation
function of a string of operators which corresponds to the domain wall states.
By identifying the correlation function of another string of operators with
the cylinder partition function of the 6-vertex model, we find that they
exhibit different behaviors in the long time limit. It would be interesting to
compute the correlation functions of more types of string operators and
investigate universal behaviors for correlation functions quantum integrable
circuits.

Another interesting problem is computing the correlators for the case where
$q$ is a root of unity, _i.e._ $q=\exp(\mathrm{i}\pi\ell_{1}/\ell_{2})$. It is
well-known that at this value, the BAE is much more subtle than the generic
$q$ case. A naive numerical solution of the BAE will lead to infinitely many
solutions due to the existence of certain bound states called Fabricius-McCoy
strings [45, 46]. The rational $Q$-system at root of unity was constructed
very recently and we now have a procedure to find all physical solutions [29].
Note that this is not only interesting from the perspective of mathematical
physics, physical systems at these values of parameters can exhibit novel
features. One example is the spin transport property at root of unity (see
_e.g_ [47] and references therein). It would be very interesting to
investigate the quantum circuit model at these values.

Finally, it would be interesting to consider other integrable quantum circuit
models such as the chiral Hubbard model in [18] with the techniques developed
in this paper. Such models correspond to spin chains with higher rank symmetry
and are technically more challenging. On the other hand, from the study of
integrable models it is expected that there would be new features such as spin
charge separation. It would be fascinating to see the effect between charge
and spin degrees of freedom in the quantum circuit models. Part of the
theoretical tools such as the rational $Q$-system have been developed in the
past years and it is now promising to study such models.

## Acknowledgements

We thank Janko Böhm, Jesper Jacobsen and Youjin Deng for the enlightening
discussions, and in particular thank Yongqun Xu for remotely operating the
computer cluster at USTC. The research was also supported by the advanced
computing resources provided by the supercomputing center of the USTC. YZ is
supported from the NSF of China through Grant No. 12075234, 11947301,
12047502, and 12247103. AH acknowledges the grants PNRR MUR Project
PE0000023-NQSTI and PRO3 Quantum Pathfinder.

##  Appendix A Proof of overlap formula (70)

In this appendix, we prove the formula (70). Consider the monodromy matrix
with generic inhomogeneity

$\displaystyle M_{0}(u|\\{\theta_{j}\\})=$
$\displaystyle\,R_{01}(u-\theta_{1})R_{02}(u-\theta_{2})\ldots
R_{0,2L}(u-\theta_{2L})=\left(\begin{array}[]{cc}A(u)&B(u)\\\
C(u)&D(u)\end{array}\right)_{0}\,.$ (120)

The Bethe states are generated by $|\mathbf{u}_{K}\rangle=B(u_{1})\ldots
B(u_{K})|\Omega\rangle$. Taking the spectral parameter at the inhomogeneities,
one notices that

$\displaystyle B(\theta_{2L-N+1})\ldots
B(\theta_{2L-1})B(\theta_{2L})|\Omega\rangle=|\underbrace{0\cdots
0}_{2L-N}\underbrace{1\cdots 1}_{N}\rangle\,.$ (121)

Similarly, for the dual state we have

$\displaystyle\langle\Omega|C(\theta_{1})C(\theta_{2})\ldots
C(\theta_{N})=\langle\underbrace{1\cdots 1}_{N}\underbrace{0\cdots 0}_{L-N}|$
(122)

Making the alternating choice of the inhomogeneities (44), it follows that we
have

$\displaystyle|0\cdots 0\underbrace{1\cdots 1}_{N}\rangle=$
$\displaystyle\,B(-\tfrac{u}{2})B(\tfrac{u}{2})\ldots
B(-\tfrac{u}{2})B(\tfrac{u}{2})|\Omega\rangle$ (123)
$\displaystyle\langle\underbrace{1\cdots 1}_{N}0\cdots 0|=$
$\displaystyle\,\langle\Omega|C(-\tfrac{u}{2})C(\tfrac{u}{2})\ldots
C(-\tfrac{u}{2})C(\tfrac{u}{2})\,.$

Using the shift operator $V$ defined in (53) and assuming that $N$ is even, we
have

$\displaystyle\langle\mathbf{u}_{N}|\text{DW}_{N}\rangle=$
$\displaystyle\,\langle\mathbf{u}_{N}|V^{N}|0\cdots 0\underbrace{1\cdots
1}_{N}\rangle$ (124) $\displaystyle=$
$\displaystyle\,\prod_{j=1}^{N}\left(\frac{\varphi(u_{j}-u/2)}{\varphi(u_{k}-u/2+\eta)}\frac{\varphi(u_{j}+u/2)}{\varphi(u_{k}+u/2+\eta)}\right)^{N/2}\langle\mathbf{u}_{N}|0\cdots
0\underbrace{1\cdots 1}_{N}\rangle$

where we have used (55) and (56). Now using the symmetric property of the
scalar product we have

$\displaystyle\langle\mathbf{u}_{N}|0\cdots 0\underbrace{1\cdots
1}_{N}\rangle=$
$\displaystyle\,\langle\Omega|\prod_{j=1}^{N}C(u_{j})B(-\tfrac{u}{2})B(\tfrac{u}{2})\ldots
B(-\tfrac{u}{2})B(\tfrac{u}{2})|\Omega\rangle$ (125) $\displaystyle=$
$\displaystyle\,\langle\Omega|C(-\tfrac{u}{2})C(\tfrac{u}{2})\ldots
C(-\tfrac{u}{2})C(\tfrac{u}{2})\prod_{j=1}^{N}B(u_{j})|\Omega\rangle$
$\displaystyle=$
$\displaystyle\,\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle\,.$

Plugging (125) to (124), we obtain

$\displaystyle\langle\mathbf{u}_{N}|\text{DW}_{N}\rangle=\prod_{j=1}^{N}\left(\frac{\varphi(u_{j}-u/2)}{\varphi(u_{k}-u/2+\eta)}\frac{\varphi(u_{j}+u/2)}{\varphi(u_{k}+u/2+\eta)}\right)^{N/2}\langle\text{DW}_{N}|\mathbf{u}_{N}\rangle,$
(126)

which is (70) in the main text.

##  Appendix B More on computational algebraic geometry

In this appendix, we give more details for the computational algebraic
geometry method (or commutative algebra) and its implementations. For a more
thorough introduction to the subject, we refer to the text book [48].

###  B.1 Polynomial ring and ideal

Let us consider a polynomial ring $\mathcal{A}=\mathbb{C}[x_{1},\ldots,x_{n}]$
which consists of all polynomials in variables $\\{x_{k}\\}$ with complex
coefficients. It is also possible to consider the polynomial ring over other
fields instead of $\mathbb{C}$.

An _ideal_ $\mathcal{I}$ of the ring $\mathcal{A}$ is a subset such that

  1. 1.

If $f_{1}\in\mathcal{I}$ and $f_{2}\in\mathcal{I}$, then
$f_{1}+f_{2}\in\mathcal{I}$;

  2. 2.

If $f\in\mathcal{I}$, then $fg\in\mathcal{I}$ for any $g\in\mathcal{A}$.

A polynomial ring is a _Noether_ ring, which means any ideal $\mathcal{I}$ of
$\mathcal{A}$ is _finitely generate_. This means for an ideal $\mathcal{I}$,
there exists a finite number of polynomials $f_{i}\in\mathcal{I}$ such that
any polynomial $F\in\mathcal{I}$ can be expressed as

$\displaystyle F=\sum_{i=1}^{K}g_{i}f_{i},\qquad g_{i}\in\mathcal{A}\,.$ (127)

These polynomials $\\{f_{1},\ldots,f_{K}\\}$ is called a basis for the ideal
$\mathcal{I}$ and we can write $\mathcal{I}=\langle
f_{1},\ldots,f_{K}\rangle$. The ideal is intimately related to the set of
polynomial equations $f_{1}=f_{2}=\ldots f_{K}=0$. In fact, we can see that if
this set of equations are satisfied, all elements in the ideal $\mathcal{I}$
is vanishing. The basis for a given ideal $\mathcal{I}$ is not unique and we
can make a ‘change of basis’ to a different set of polynomials

$\displaystyle\mathcal{I}=\langle f_{1},\ldots,f_{K}\rangle=\langle
g_{1},\ldots,g_{K^{\prime}}\rangle\,,$ (128)

where in general $K\neq K^{\prime}$.

###  B.2 Gröbner basis and quotient ring

Given an ideal $\mathcal{I}$, we can define the quotient ring
$V_{f}=\mathcal{A}/\mathcal{I}$ as the quotient set specified by the following
equivalence relation

$\displaystyle F\sim G\qquad\text{if and only if}\qquad F-G\in\mathcal{I}\,.$
(129)

This means for a given basis of the ideal $\mathcal{I}=\langle
f_{1},\ldots,f_{k}\rangle$ we make the following identification

$\displaystyle F\sim F+\sum_{i=1}^{K}a_{i}f_{i},\qquad a_{k}\in\mathcal{A}\,.$
(130)

Given a polynomial $F(u)$, we can perform a polynomial reduction with respect
to $f_{i}$ and find a remainder. Roughly speaking, the quotient ring is the
space of such remainders. However, for an arbitrary generating set, the
aforementioned polynomial reduction is ill-defined because the ‘remainder’ is
not unique. Therefore, in order to make the reduction in a polynomial ring
possible, one needs a “convenient” set of basis such that the polynomial
reduction leads to a well-defined and unique remainder. Such convenient basis
indeed exists and are called Gröbner basis of the ideal.

To define a Gröbner basis, we need to first define monomial orders in a
polynomial ring. A monomial orders $\prec$ is a total order for all monomials
in $\mathcal{A}$ such that

  * •

If $u\prec v$ the for any monomial $w$, $uw\prec vw$.

  * •

If $u$ is a non-constant monomial, then $1\prec u$.

In most of our calculations, we choose the _degrevlex_
(DegreeReverseLexicographic) ordering. After defining a monomial order
$\prec$, for any polynomial $F\in\mathcal{A}$, there is a unique _leading
term_ , denoted by $\text{LT}(f)$ which is the highest monomial of $F$ in the
order $\prec$. Now we can give the formal definition of Gröbner basis.

#### Gröbner basis

A Gröbner basis $G(\mathcal{I})$ of an ideal $\mathcal{I}$ with respect to a
monomial order $\prec$ is a generating set of $\mathcal{I}$ such that for any
$F\in\mathcal{I}$,

$\displaystyle\exists g_{i}\in
G(\mathcal{I}),\qquad\text{LT}(g_{i})|\text{LT}(F)\,$ (131)

where $a|b$ means a monomial $b$ is divisible by another monomial $a$. Given a
monomial order $\prec$ and any generating set $f_{i}$ of the ideal
$\mathcal{I}$, there are standard algorithms such as the Buchberger algorithm
or the F4/F5 algorithm to compute the Gröbner basis. Such algorithms have been
implemented in many commonly used softwares such as Mathematica or more
specialized packages such as Singular. For example, we can use the Mathematica
function GroebnerBasis to compute the Gröbner basis of a given ideal.

The property (131) ensures that the polynomial division of any polynomial
$F\in\mathcal{A}$ towards an ideal $\mathcal{I}$ in the order $\prec$ is well-
defined

$\displaystyle F=\sum a_{i}g_{i}+r$ (132)

where $g_{i}$’s are the elements of the Gröbner basis and $t$ is called the
remainder, which contains monomials _not_ divisible by any $\text{LT}(g_{i})$.
The key point is that, given the monomial order $\prec$, the remainder $r$ for
$F$ is unique.

#### Quotient ring basis

With Gröbner basis, the polynomial reduction now can be performed. The
quotient ring $V_{f}=\mathcal{A}/\mathcal{I}$ can now be identified with the
space of ‘remainders’. A canonical basis for the quotient ring can be chosen
as the monomials which are not divisible by any elements in
$\text{LT}(G(\mathcal{I}))$. Let us denote this basis by
$B=\\{b_{1},b_{2},\ldots,b_{N_{s}}\\}$. One important result from algebraic
geometry is that the number of solutions of the equation
$f_{1}=f_{2}\ldots=f_{K}=0$ equals the dimension of the quotient ring $\dim
V_{f}=N_{s}$. This statement provides a valuable method for determining the
number of solutions. In practice, we can use the lattice algorithm to list
these monomials.

###  B.3 Companion matrix

Let $B=\\{b_{1},\ldots,b_{N_{s}}\\}$ be the canonical basis of $V_{f}$, for
any $F\in\mathcal{A}$, we can define

$Fb_{i}\sim\sum_{j=1}^{N_{s}}b_{j}(\mathbf{M}_{F})_{ji}$ (133)

where $\sim$ stands for the procedure of performing polynomial reduction with
respect to Gröbner basis and finding the remainder. $\mathbf{M}_{F}$ is an
$N_{s}\times N_{s}$ matrix which is called the companion matrix of $F$. Let us
denote the solutions of the equations $f_{1}=f_{2}=\ldots=f_{K}=0$ by
$\mathcal{Z}(I)$. We have the following important result

$\displaystyle\sum_{p\in\mathcal{Z}(I)}f(p)=\text{tr}(\mathbf{M}_{F})$ (134)

This is the key formula for our purpose. Note that to calculate the right hand
side, we only need to compute $G(\mathcal{I})$ which is purely algebraic and
analytic. The companion matrices have the following properties

$\displaystyle\mathbf{M}_{F_{1}+F_{2}}=\mathbf{M}_{F_{1}}+\mathbf{M}_{F_{2}},\qquad\mathbf{M}_{F_{1}F_{2}}=\mathbf{M}_{F_{1}}\mathbf{M}_{F_{2}}=\mathbf{M}_{F_{2}}\mathbf{M}_{F_{1}}$
(135)

and also

$\displaystyle\mathbf{M}_{F_{1}/F_{2}}=\mathbf{M}_{F_{1}}\cdot\mathbf{M}_{F_{2}}^{-1}\,.$
(136)

This last property allows us to compute the companion matrix of a rational
function $f/g$ as

$\displaystyle\sum_{p\in\mathcal{Z}(I)}\frac{f(p)}{g(p)}=\text{tr}(M_{f}M_{g}^{-1})\,,$
(137)

as long as $g(p)\not=0$, $\forall p\in\mathcal{Z}(\mathcal{I})$.

###  B.4 Efficient implementations

In the computational algebraic geometry approach, the majority of the
computational time is spent on the computation of the Gröbner basis $G(I)$ and
the companion matrix [30]. In this work, the former computation is powered by
the Buchberger algorithm implemented in Singular [49] and also the new
parallel Buchberger algorithm implemented with Singular/GPI-space framework
[50] [51]. The latter computation is powered by the multivariate polynomial
division function in Singular. In practice, the coefficients of the equations
under consideration are not rational or complex numbers, but rational
functions of some physical parameters under certain algebraic extension. The
Gröbner basis computation with parameters is possible but inefficient. We
carry out the Gröbner basis computation with integer-valued parameters for
better performance. The same computation with different values of parameters
can be performed in parallelization. The final physical results will come from
the interpolation of the different parameters.

In the context of this work, the results are presented in the form of rational
functions. For XXX model, the results are univariate rational functions. For
XXZ model, they are multivariate (bivariate) rational functions. To finish
Gröbner bases computation with higher efficiency, one can perform integer-
value seeding over these variables on the quotient ring representation, which
will give a large amount of simpler CAG results with fixed values of
variables.

Then to get the full analytic result in the form of rational functions, one
can perform Thiele’s rational interpolation formula for dense rational
reconstruction [52].

This algorithm is superiorly efficient for univariate interpolation of
rational function and sensitive to the number of numerical points as probes in
the computation. This process is finished with private implementation in
Singular and Mathematica.

##  Appendix C Rational $Q$-system

In this appendix, we introduce the rational $Q$-system. We present the
construction for XXZ spin chain. The rational $Q$-system for XXX is similar.

#### Young diagram

We start with a two-row Young diagram with boxes $(M,L-M)$ where $M$ is the
number of magnons and $L$ is the length of the spin chain. We consider the
regime $M\leq L/2$. At each node of the Young diagram, we associate with a
$Q$-function denoted by $Q_{a,s}(x)$ where $a$ and $s$ are the lattice
coordinates. For XXZ spin chain, we consider the $Q$-functions as a function
of the multiplicative variable $x=e^{2u}$.

####  $QQ$-relations

The $Q$-functions around each box whose lower-left corner coordinate is
$(a,s)$ are related by the $QQ$-relations

$\displaystyle
Q_{a+1,s}Q_{a,s+1}=Q_{a+1,s+1}^{+}Q_{a,s}^{-}-Q_{a+1,s+1}^{-}Q_{a,s}^{+}$
(138)

where

$\displaystyle Q^{\pm}_{a,s}(x)=Q_{a,s}(xq^{\pm 1})\,.$ (139)

Among all the $Q$-functions, the $Q_{1,0}(x)$ is special because its zeros are
the Bethe roots. We stress that the $QQ$-relation is defined up to
proportionality. So the global normalizations of the $Q$-functions are not
important.

The $QQ$-relation clearly states that the $Q$-functions at different nodes are
not independent. As a result, we can work out all the $Q$-functions starting
with a few fixed $Q$-functions. These fixed ones are the boundary conditions.

#### Boundary conditions

We fix the $Q$-functions at the upper and left boundary. The $Q$-functions at
the upper boundary are fixed to be a constant $Q_{2,s}(x)=1$. The rest two
$Q$-functions at the left boundary is given by

$\displaystyle Q_{0,0}(x)=(x^{1/2}-x^{-1/2})^{L},\qquad
Q_{1,0}(x)=\prod_{j=1}^{M}\left((x/x_{j})^{1/2}-(x_{j}/x)^{1/2}\right)$ (140)

where $\\{x_{j}\\}$ are the Bethe roots. Since the global normalization is not
important, we have the following parametrization

$\displaystyle
Q_{1,0}(x)\propto\widetilde{Q}_{1,0}(x)=\prod_{j=1}^{M}(x-x_{j})=x^{M}+\sum_{k=1}^{M-1}c_{k}x^{k}\,.$
(141)

#### Zero remainder condition

To find the Bethe roots, we need to find $Q_{1,0}$, or equivalently the
coefficients $\\{c_{k}\\}$. This is achieved by the zero remainder conditions.
By fixing the upper and boundary conditions, we can work out the rest of the
$Q$-functions by the $QQ$-relations. However, for a generic choice of
$\\{c_{k}\\}$, the resulting $Q$-functions are not Laurent polynomials in $x$
(up to an unimportant global factor). By _requiring_ all the $Q$-functions to
be Laurent polynomials in $x$, we obtain further constraints on $\\{c_{j}\\}$,
which takes the form of a set of algebraic equations. These equations are
called the _zero remainder conditions_.

One important point about the zero remainder condition is that, in fact we
only need to require the $Q$-functions of the first two columns to be Laurent
polynomials. More precisely, we only need to make sure that $Q_{0,1}$ and
$Q_{0,2}$ are polynomials. It can be proven that this is sufficient to
guarantee that the rest of the $Q$-functions are Laurent polynomials.

##  Appendix D Symmetrization

In this appendix, we present more details for the change of variables from
$\\{x_{k}\\}$ to $\\{c_{k}\\}$. In general, we have an expression of the form
$N(\\{x_{j}\\})/D(\\{x_{j}\\})$ where $N(\\{x_{j}\\})$ and $D(\\{x_{j}\\})$
are polynomials in $\\{x_{j}\\}$. We first establish the relation between
$\\{x_{j}\\}$ and $\\{c_{j}\\}$ by identifying

$\displaystyle\prod_{j=1}^{M}(x-x_{j})=x^{M}+\sum_{k=1}^{M-1}c_{k}\,x^{k}\,.$
(142)

These lead to

$\displaystyle c_{M-1}=$ $\displaystyle\,-(x_{1}+x_{2}+\ldots+x_{M})\,,$ (143)
$\displaystyle c_{M-2}=$
$\displaystyle\,x_{1}x_{2}+x_{1}x_{3}+\ldots+x_{M-1}x_{M}\,,$
$\displaystyle\cdots$ $\displaystyle c_{0}=$
$\displaystyle\,(-1)^{M}x_{1}x_{2}\ldots x_{M}\,.$

We then introduce another equation

$\displaystyle
w=\frac{N(\\{x_{j}\\})}{D(\\{x_{j}\\})}\quad\Leftrightarrow\quad
N(\\{x_{j}\\})=w\,D(\\{x_{j}\\})\,.$ (144)

Combing the two sets of equations (143) and (144), we obtain a system of
algebraic equation involving
$\\{w,x_{1},\ldots,x_{M},c_{0},\ldots,c_{M-1}\\}$. We can then compute the
Gröbner basis of these set of algebraic equations. The key observation is
that, among the Gröbner basis, there is always a polynomial that involve only
$\\{w,c_{0},\ldots,c_{M-1}\\}$ and the dependence on $w$ is linear. Let us
denote this polynomial by $S(w,c_{0},\ldots,c_{M-1})$. By solving the equation

$\displaystyle S(w,c_{0},\ldots,c_{M-1})=0\,.$ (145)

for $w$. Since the dependence on $w$ is linear, it can be solve readily and
obtain $w=\widetilde{S}(\\{c_{k}\\})$. The quantity
$\widetilde{S}(\\{c_{k}\\})$ is the result of symmetrization of
$N(\\{x_{j}\\})/D(\\{x_{j}\\})$, _i.e._

$\displaystyle\widetilde{S}(\\{c_{k}\\})=\frac{N(\\{x_{j}\\})}{D(\\{x_{j}\\})}\,.$
(146)

##  Appendix E More analytic results in massive regime

In this appendix, we present analytic results for $L=6,\ldots,10$ and $M=2$
and $L=8$, $M=4$ in Fourier space. We fix $\xi=9$, $q=2$ as in the main text.

#### Result for $L=6,M=2$

For $L=6$, $M=2$, we have

$\displaystyle\hat{P}_{6}^{(2)}(a)=$
$\displaystyle\,-5341150064515340273121\,a^{6}-169734594354674880220074\,a^{5}\,$
$\displaystyle\,-67729782571339869959686015\,a^{4}+1208144728170975147670407220\,a^{3}$
$\displaystyle\,+1370616640991778706771271985\,a^{2}+113499929340996879020991717526\,a$
$\displaystyle\,-1289526489840334121941717217521\,,$
$\displaystyle\hat{Q}_{6}^{(2)}(a)=$
$\displaystyle\,43439888521963583647921\,a^{7}-1100559700418854159130647\,a^{6}$
$\displaystyle\,-1768438805199882583288059\,a^{5}-151648502182721393666677235\,a^{4}$
$\displaystyle\,+2397197990231124981189565235\,a^{3}+2236019003611924904658601659\,a^{2}$
$\displaystyle\,+122510509439470012808294498647\,a$
$\displaystyle-1289526489840334121941717217521.$

#### Result for $L=7,M=2$

For $L=7$, $M=2$, we have

$\displaystyle\hat{P}_{7}^{(2)}(a)$
$\displaystyle=\,-31667678732511452479334409\,a^{7}-972253548728081689402947537\,a^{6}$
$\displaystyle\,-13284014090995839018704234589\,a^{5}-2658877056547083933342929909685\,a^{4}$
$\displaystyle\,-7056089646281879410086216120315\,a^{3}-132402145030090476007173235599411\,a^{2}$
$\displaystyle\,-12309244542150474755821891018602463\,a$
$\displaystyle+134535009158552218608057415586748409$
$\displaystyle\hat{Q}_{7}^{(2)}a)$
$\displaystyle=\,257555099046722087448523609\,a^{8}-6802579045483603414528506472\,a^{7}$
$\displaystyle\,-9072737719011528891925528548\,a^{6}-192957938169692945243642474904\,a^{5}$
$\displaystyle\,-5305531949916488629095660263370\,a^{4}-10904271570463836869480995222104\,a^{3}$
$\displaystyle\,-219419179113019027067093413101348\,a^{2}$
$\displaystyle-13249309353244078330717402870175272\,a$
$\displaystyle\,+134535009158552218608057415586748409\,.$

#### Result for $L=8,M=2$

For $L=8$, $M=2$, we have

$\displaystyle\hat{P}_{8}^{(2)}(a)=$
$\displaystyle\,-187757667205060401749973710961\,a^{8}$
$\displaystyle-5562295426349338066893801548712\,a^{7}$
$\displaystyle\,-68460063735650825294029417860108\,a^{6}$
$\displaystyle-43586334407867250947331928333142584\,a^{5}$
$\displaystyle\,+777575727956544615801669787675266730\,a^{4}$
$\displaystyle\,+626223846698148298352660070921579816\,a^{3}$
$\displaystyle\,+12538451422031669290868599967103417492\,a^{2}$
$\displaystyle\,+1333027799413296620840891103020100760888\,a$
$\displaystyle\,-14035902970502594415160022110749874762561$
$\displaystyle\hat{Q}_{8}^{(2)}(a)=$
$\displaystyle\,1527044182248015256482296477761\,a^{9}$
$\displaystyle-41976962049572871860538871397449\,a^{8}$
$\displaystyle\,-43647594567560624615959514672604a^{7}$
$\displaystyle-950080126167093004906431342935124\,a^{6}$
$\displaystyle\,-97727582196713429284273985905005314\,a^{5}$
$\displaystyle\,+1543505891415224072829997935348218114\,a^{4}$
$\displaystyle\,+973825436919024486110553161407549524\,a^{3}$
$\displaystyle\,+21275744879111919409721635343616176604\,a^{2}$
$\displaystyle\,+1431103821089881188206164958982840351049\,a$
$\displaystyle\,-14035902970502594415160022110749874762561\,.$

#### Result for $L=9,M=2$

For $L=9$, $M=2$, we have

$\displaystyle\hat{P}_{9}^{(2)}(a)=$
$\displaystyle\,-1113215208858803121975594132287769\,a^{9}$
$\displaystyle\,-31780030304816697318295618475616879\,a^{8}$
$\displaystyle\,-346116949924399518400809244763290884\,a^{7}$
$\displaystyle\,-3421179821316113328851446955867198604\,a^{6}$
$\displaystyle\,-2172371663730681721093427279453411135694\,a^{5}$
$\displaystyle\,-2857540521689186032968321538208467721106\,a^{4}$
$\displaystyle\,-53168876840828321946403512571977219887796\,a^{3}$
$\displaystyle\,-1157398391847265559110526743401755305933116\,a^{2}$
$\displaystyle\,-144166447014633183154420633953273196460154721\,a$
$\displaystyle\,+1464351721009565172739229946792423684103226569\,,$
$\displaystyle\hat{Q}_{9}^{(2)}(a)=$
$\displaystyle\,9053844956548482455683535816644969\,a^{10}$
$\displaystyle\,-258631475892209138863609353351968890\,a^{9}$
$\displaystyle\,-188139026338731523604678897771121195\,a^{8}$
$\displaystyle\,-4335214914789008029067449828040746680\,a^{7}$
$\displaystyle\,-108767060827438760278145243348349756510\,a^{6}$
$\displaystyle\,-4342360265037034821500357626239933445788\,a^{5}$
$\displaystyle\,-3833731890336342622524151459421149169310\,a^{4}$
$\displaystyle\,-83936810522729090025396186967817696823480\,a^{3}$
$\displaystyle\,-2033364160886612043270639024726458464777995\,a^{2}$
$\displaystyle\,-154398620280129574483072290072009855162061690\,a$
$\displaystyle\,+1464351721009565172739229946792423684103226569$

#### Result for $L=10,M=2$

For $L=10$, $M=2$, we have

$\displaystyle\hat{P}_{10}^{(2)}(a)=$
$\displaystyle\,-6600252973323843710193297610334182401\,a^{10}$
$\displaystyle\,-181316000177945635411970895396128517990\,a^{9}$
$\displaystyle\,-1705329740200620956534465254659507286845\,a^{8}$
$\displaystyle\,-10754332693422987296905079249680362609480\,a^{7}$
$\displaystyle\,-27885978656919283770966931750762365032046610\,a^{6}$
$\displaystyle\,+497564943588073965118725621401742897822373852\,a^{5}$
$\displaystyle\,+200286748319860161399864219482404031581811790\,a^{4}$
$\displaystyle\,+4211824928912551127213750961543388449686740920\,a^{3}$
$\displaystyle\,+103177140640554805178621133752236693674989553555\,a^{2}$
$\displaystyle\,+15572087776093627470623102170300915263222808880410\,a$
$\displaystyle\,-152774350701206924906711121118906770538805524717201\,.$
$\displaystyle\hat{Q}_{10}^{(2)}(a)=$
$\displaystyle\,53680246747375952479747683856888021201a^{11}$
$\displaystyle\,-1591234173145736771643410483719395105211\,a^{10}$
$\displaystyle\,-634353406652137645111028847889959009945\,a^{9}$
$\displaystyle\,-17199644191884346788964056757783189794165\,a^{8}$
$\displaystyle\,-477986240728410409714941746696418174059670\,a^{7}$
$\displaystyle\,-62641723481738385823504462469886633388210862\,a^{6}$
$\displaystyle\,+988211694405938127147292885147793827148962862\,a^{5}$
$\displaystyle\,+252685292016123719737141508261795875216379670\,a^{4}$
$\displaystyle\,+6861229955645771169924139291085012167575594165\,a^{3}$
$\displaystyle\,+190852983845750775085500945828236892588740889945\,a^{2}$
$\displaystyle\,+16639600180709600481550000801512592128934061049211\,a$
$\displaystyle\,-152774350701206924906711121118906770538805524717201$

We also present the result for $L=4$, $M=4$. As we can see, the result is
considerably more complicated.

#### Result for $L=4,M=4$

$\displaystyle\hat{P}_{4}^{(4)}(a)=\,-8156172132531761716980920529095901259254650042562081a^{11}$
$\displaystyle\,-20998071783853815347919633492950500857096693679074341509\,a^{10}$
$\displaystyle\,+9805779563942996214406401795772195686233369021819555215945\,a^{9}$
$\displaystyle\,-486409425272801521035561455424640295699958396718871092832235\,a^{8}$
$\displaystyle\,-19398767636794702594635762313443606166288950348747029297528330\,a^{7}$
$\displaystyle\,+1379499453838511260832184210395615109863035010775248140015259822\,a^{6}$
$\displaystyle\,+1049908096848413039112303771545015356006611044880065645611673778\,a^{5}$
$\displaystyle\,-1067502230719708323373110373611855273607688310911791538353875523670\,a^{4}$
$\displaystyle\,+13757920624322051646768038590766100415496552222863074579585797842235\,a^{3}$
$\displaystyle\,+140206788726046972218515610981732164173910159246224337801102059886055\,a^{2}$
$\displaystyle\,-362405203824894683955033973028680510437736739186394290704287512112129\,a$
$\displaystyle\,+1809966868929172933994874978249710316608864595987799100807916546403128$
$\displaystyle\hat{Q}_{4}^{(4)}(a)=\,66334629044966824718765264495212069363677393830384881\,a^{12}$
$\displaystyle\,-27643127263788309609655641814912881709982050118818188172\,a^{11}$
$\displaystyle\,+1360610363855872619271697002210675787787122935496392344546\,a^{10}$
$\displaystyle\,+73845237226000279421525338791921188138471306476200512550180\,a^{9}$
$\displaystyle\,-4780979055224263220320109492155092475298551712325064333203905\,a^{8}$
$\displaystyle\,-37562545606540651114101312071470418979689832312217881666249752\,a^{7}$
$\displaystyle\,+5214641940831811898196292622170202201005928496707287248866193244\,a^{6}$
$\displaystyle\,-42518415110005549339013514206450781815189905693149656045060534552\,a^{5}$
$\displaystyle\,-1658663539373636328709728248947847374053887035751846492027459307905\,a^{4}$
$\displaystyle\,+26300083629274148442984818223919740731189116203000185638476277686180\,a^{3}$
$\displaystyle\,+95689865821052118423630225085462686944116781710236380669755943826146\,a^{2}$
$\displaystyle\,-375052366621905745777559482322779359471158532652018475202259994953217\,a$
$\displaystyle\,+1809966868929172933994874978249710316608864595987799100807916546403128\,.$

##  Appendix F More analytic results in massless regime

In this appendix, we present analytic results for $L=4,\ldots,8$ and $M=2$ and
$L=8$, $M=4$ in Fourier space. We fix $\xi=9$,
$q=\tfrac{3}{5}+\tfrac{4\mathrm{i}}{5}$ as in the main text.

#### Result for $L=4,M=2$

$\displaystyle\hat{P}_{4}^{(2)}(a)=$
$\displaystyle\,(16190945280+23440752896\,\mathrm{i})a^{4}+(3909369240-18626610017\,\mathrm{i})a^{3}$
$\displaystyle\,+(18961677720+16502505825\,\mathrm{i})a^{2}-(32861123280+17495261183\,\mathrm{i})a$
$\displaystyle\,+(12788123040+51977892479\,\mathrm{i})\,.$
$\displaystyle\hat{Q}_{4}^{(2)}(a)=$
$\displaystyle\,(12788123040-51977892479\,\mathrm{i})a^{5}-(16670178000-40936014079\,\mathrm{i})a^{4}$
$\displaystyle\,+(22871046960-35129115842\,\mathrm{i})a^{3}+(22871046960+35129115842\,\mathrm{i})a^{2}$
$\displaystyle\,-(16670178000+40936014079\,\mathrm{i})a+(12788123040+51977892479\,\mathrm{i})\,.$

#### Result for $L=5,M=2$

$\displaystyle\hat{P}_{5}^{(2)}(a)=$
$\displaystyle\,-(13306340474624+3273759498240\,\mathrm{i})a^{5}$
$\displaystyle\,+(3863467459072+13489014988800\,\mathrm{i})a^{4}$
$\displaystyle\,+(7773196913152-3163367669760\,\mathrm{i})a^{3}$
$\displaystyle\,-(6745452149248+12433548165120\,\mathrm{i})a^{2}$
$\displaystyle\,+(10796522125247+19922035736520\,\mathrm{i})a$
$\displaystyle\,+(11977223406401-22791452542200\,\mathrm{i})\,.$
$\displaystyle\hat{Q}_{5}^{(2)}(a)=$
$\displaystyle\,(11977223406401+22791452542200\,\mathrm{i})a^{6}$
$\displaystyle\,-(2509818349377+23195795234760\,\mathrm{i})a^{5}$
$\displaystyle\,-(2881984690176-25922563153920\,\mathrm{i})a^{4}$
$\displaystyle\,+15546393826304a^{3}$
$\displaystyle\,-(2881984690176+25922563153920\,\mathrm{i})a^{2}$
$\displaystyle\,-(2509818349377-23195795234760\,\mathrm{i})a$
$\displaystyle\,+(11977223406401-22791452542200\,\mathrm{i})\,.$

#### Result for $L=6,M=2$

$\displaystyle\hat{P}_{6}^{(2)}=$
$\displaystyle\,(3066169192038656+5834611850803200\,\mathrm{i})a^{6}$
$\displaystyle\,+(4787760828350976-6497951173079040\,\mathrm{i})a^{5}$
$\displaystyle\,-(2879859037024545+105889836591000\,\mathrm{i})a^{4}$
$\displaystyle\,+(2599750449819425-2871676258248600\,\mathrm{i})a^{3}$
$\displaystyle\,-(7491925945812480-133968233779200\,\mathrm{i})a^{2}$
$\displaystyle\,+(13197453475178049+1833272092202400\,\mathrm{i})a$
$\displaystyle\,-(4384188648550081+11582273787266160\,\mathrm{i})\,.$
$\displaystyle\hat{Q}_{6}^{(2)}=$
$\displaystyle\,(4384188648550081-11582273787266160\,\mathrm{i})a^{7}$
$\displaystyle\,-(10131284283139393-7667883943005600\,\mathrm{i})a^{6}$
$\displaystyle\,+(12279686774163456-6363982939299840\,\mathrm{i})a^{5}$
$\displaystyle\,-(5479609486843970+2977566094839600\,\mathrm{i})a^{4}$
$\displaystyle\,+(5479609486843970-2977566094839600\,\mathrm{i})a^{3}$
$\displaystyle\,-(12279686774163456+6363982939299840\,\mathrm{i})a^{2}$
$\displaystyle\,+(10131284283139393+7667883943005600\,\mathrm{i})a$
$\displaystyle\,-(4384188648550081+11582273787266160\,\mathrm{i})\,.$

#### Result for $L=7,M=2$

$\displaystyle\hat{P}_{7}^{(2)}=$
$\displaystyle\,-(2965062089540136960+1122352294028820736\,\mathrm{i})a^{7}$
$\displaystyle\,+(186664534285271040+4527404219884142080\,\mathrm{i})a^{6}$
$\displaystyle\,+(1883104092849070080-1055138737458676224\,\mathrm{i})a^{5}$
$\displaystyle\,-(287336924931686400+531319906398277120\,\mathrm{i})a^{4}$
$\displaystyle\,-(2502921301220966400-653553357493434880\,\mathrm{i})a^{3}$
$\displaystyle\,+(3740222195766927360+1789593583871064576\,\mathrm{i})a^{2}$
$\displaystyle\,-(5298184152707352840+5464802435030160895\,\mathrm{i})a$
$\displaystyle\,+(-2116437424659875880+5568174742303293439\,\mathrm{i})\,.$
$\displaystyle\hat{Q}_{7}^{(2)}=$
$\displaystyle\,(2116437424659875880+5568174742303293439\,\mathrm{i})a^{8}$
$\displaystyle\,+(2333122063167215880-6587154729058981631\,\mathrm{i})a^{7}$
$\displaystyle\,-(3553557661481656320-6316997803755206656\,\mathrm{i})a^{6}$
$\displaystyle\,+(4386025394070036480-401585379965241344\,\mathrm{i})a^{5}$
$\displaystyle\,-(1062639812796554240\,\mathrm{i})a^{4}$
$\displaystyle\,-(4386025394070036480+401585379965241344\,\mathrm{i})a^{3}$
$\displaystyle\,+(3553557661481656320+6316997803755206656\,\mathrm{i})a^{2}$
$\displaystyle\,-(2333122063167215880+6587154729058981631\,\mathrm{i})a$
$\displaystyle\,(-2116437424659875880+5568174742303293439\,\mathrm{i})\,.$

#### Result for $L=8,M=2$

$\displaystyle\hat{P}_{8}^{(2)}=$
$\displaystyle\,(541807980712928225280+1425452734029643120384\,\mathrm{i})a^{8}$
$\displaystyle\,+(1860089695679579811840-1663719405239146709504\,\mathrm{i})a^{7}$
$\displaystyle\,-(1306086559669141401600+583178930942280794624\,\mathrm{i})a^{6}$
$\displaystyle\,+(39299130340212338280+290838133169585319649\,\mathrm{i})a^{5}$
$\displaystyle\,-(283996239369855701400+243181411869805030625\,\mathrm{i})a^{4}$
$\displaystyle\,-(555465568612198625280-1485219409999454373376\,\mathrm{i})a^{3}$
$\displaystyle\,+(1926211061318640168960-1440065866208923122176\,\mathrm{i})a^{2}$
$\displaystyle\,-(4232671666695290048400-381810663241616919679\,\mathrm{i})a$
$\displaystyle\,(1329399368762685232320+2538165215672305923841\,\mathrm{i})\,.$
$\displaystyle\hat{Q}_{8}^{(2)}=$
$\displaystyle\,(1329399368762685232320-2538165215672305923841\,\mathrm{i})a^{9}$
$\displaystyle\,-(3690863685982361823120-1043642070788026200705\,\mathrm{i})a^{8}$
$\displaystyle\,+(3786300756998219980800-223653539030223587328\,\mathrm{i})a^{7}$
$\displaystyle\,-(1861552128281340026880+2068398340941735168000\,\mathrm{i})a^{6}$
$\displaystyle\,-(244697109029643363120-534019545039390350274\,\mathrm{i})a^{5}$
$\displaystyle\,-(244697109029643363120+534019545039390350274\,\mathrm{i})a^{4}$
$\displaystyle\,-(1861552128281340026880-2068398340941735168000\,\mathrm{i})a^{3}$
$\displaystyle\,+(3786300756998219980800+223653539030223587328\,\mathrm{i})a^{2}$
$\displaystyle\,-(3690863685982361823120+1043642070788026200705\,\mathrm{i})a$
$\displaystyle\,+(1329399368762685232320+2538165215672305923841\,\mathrm{i})\,.$

#### Result for $L=4,M=4$

$\displaystyle\hat{P}_{4}^{(4)}$
$\displaystyle=\,-(25914677165518638552937393664573696+18885372235073878041827729279659653120\,\mathrm{i})a^{11}$
$\displaystyle\,+(24630256608315177583960222249037046721+14791618858230018600341357060533136880\,\mathrm{i})a^{10}$
$\displaystyle\,-(35277224091134454445608489237612217220-10852527023685678060739087229387853360\,\mathrm{i})a^{9}$
$\displaystyle\,+(15889997309911171224836659068505377795-29641482629099231105353029850095050160\,\mathrm{i})a^{8}$
$\displaystyle\,-(2354652261433981164630005385296711160-26442493434928705349377006828321757760\,\mathrm{i})a^{7}$
$\displaystyle\,-(24084748945745632708905819811178781390+238832243927563314220178586885848160\,\mathrm{i})a^{6}$
$\displaystyle\,-(7778950403883894847540608617741706240+25918642395117628948049529404717382240\,\mathrm{i})a^{5}$
$\displaystyle\,+(49077616940393180206610546519782250190+29836741776308322341792887609402796640\,\mathrm{i})a^{4}$
$\displaystyle\,-(42996910595483003653184616761086396680-14211989812263351882467726396115543360\,\mathrm{i})a^{3}$
$\displaystyle\,+(45532442407233332355320615831884779005-48519927643957199603125363113009532560\,\mathrm{i})a^{2}$
$\displaystyle\,+(3288660014426965141749600884425604996+41361215880861563790520315012683080880\,\mathrm{i})a$
$\displaystyle\,(-26525262588323420505310919969134672321-23569386823313100801648840241916702640\,\mathrm{i})\,.$
$\displaystyle\hat{Q}_{4}^{(4)}$
$\displaystyle=\,-(26525262588323420505310919969134672321-23569386823313100801648840241916702640\,\mathrm{i})a^{12}$
$\displaystyle\,+(3262745337261446503196663490761031300-60246588115935441832348044292342734000\,\mathrm{i})a^{11}$
$\displaystyle\,+(70162699015548509939280838080921825726+63311546502187218203466720173542669440\,\mathrm{i})a^{10}$
$\displaystyle\,-(78274134686617458098793105998698613900+3359462788577673821728639166727690000\,\mathrm{i})a^{9}$
$\displaystyle\,+(64967614250304351431447205588287627985-59478224405407553447145917459497846800\,\mathrm{i})a^{8}$
$\displaystyle\,-(10133602665317876012170614003038417400-52361135830046334297426536233039140000\,\mathrm{i})a^{7}$
$\displaystyle\,-48169497891491265417811639622357562780a^{6}$
$\displaystyle\,-(10133602665317876012170614003038417400+52361135830046334297426536233039140000\,\mathrm{i})a^{5}$
$\displaystyle\,+(64967614250304351431447205588287627985+59478224405407553447145917459497846800\,\mathrm{i})a^{4}$
$\displaystyle\,-(78274134686617458098793105998698613900-3359462788577673821728639166727690000\,\mathrm{i})a^{3}$
$\displaystyle\,+(70162699015548509939280838080921825726-63311546502187218203466720173542669440\,\mathrm{i})a^{2}$
$\displaystyle\,+(3262745337261446503196663490761031300+60246588115935441832348044292342734000\,\mathrm{i})a$
$\displaystyle\,+(-26525262588323420505310919969134672321-23569386823313100801648840241916702640\,\mathrm{i})\,.$

## References

  * [1] B. Sutherland,  _Beautiful Models_ ,  World Scientific Publishing Company (2004).
  * [2] V. Korepin, N. Bogoliubov and A. Izergin,  _Quantum inverse scattering method and correlation functions_ ,  Cambridge University Press (1993).
  * [3] A. H. van Amerongen, J. J. P. van Es, P. Wicke, K. V. Kheruntsyan and N. J. van Druten,  _Yang-Yang Thermodynamics on an Atom Chip_ ,  Phys. Rev. Lett. 100(9), 090402 (2008),  10.1103/PhysRevLett.100.090402,  0709.1899.
  * [4] B. Lake, D. A. Tennant, J.-S. Caux, T. Barthel, U. Schollwöck, S. E. Nagler and C. D. Frost,  _Multispinon Continua at Zero and Finite Temperature in a Near-Ideal Heisenberg Chain_ ,  Phys. Rev. Lett. 111(13), 137205 (2013),  10.1103/PhysRevLett.111.137205,  1307.4071.
  * [5] R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl and K. Kiefer,  _Quantum Criticality in an Ising Chain: Experimental Evidence for Emergent E 8 Symmetry_,  Science 327, 177 (2010),  10.1126/science.1180085,  1103.3694.
  * [6] J. Wu, M. Kormos and Q. Si,  _Finite-temperature spin dynamics in a perturbed quantum critical ising chain with an ${E}_{8}$ symmetry_,  Phys. Rev. Lett. 113, 247201 (2014),  10.1103/PhysRevLett.113.247201.
  * [7] Z. Zhang _et al._ ,  _Observation of E8 Particles in an Ising Chain Antiferromagnet_ ,  Phys. Rev. B 101(22), 220411 (2020),  10.1103/PhysRevB.101.220411,  2005.13772.
  * [8] 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 _et al._ ,  _${E}_{8}$ Spectra of Quasi-One-Dimensional Antiferromagnet ${\mathrm{BaCo}}_{2}{\mathrm{V}}_{2}{\mathrm{O}}_{8}$ under Transverse Field_,  Phys. Rev. Lett. 127, 077201 (2021),  10.1103/PhysRevLett.127.077201.
  * [9] J. Yang, X. Wang and J. Wu,  _Magnetic excitations in the one-dimensional Heisenberg-Ising model with external fields and their experimental realizations_ ,  Journal of Physics A Mathematical General 56(1), 013001 (2023),  10.1088/1751-8121/acad48,  2301.05823.
  * [10] T. Kinoshita, T. Wenger and D. S. Weiss,  _A quantum Newton’s cradle_ ,  Nature 440, 900 (2006),  10.1038/nature04693.
  * [11] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer and J. Schmiedmayer,  _Experimental observation of a generalized Gibbs ensemble_ ,  Science 348, 207 (2015),  10.1126/science.1257026,  1411.7185.
  * [12] M. Schemmer, I. Bouchoule, B. Doyon and J. Dubail,  _Generalized HydroDynamics on an Atom Chip_ ,  Phys. Rev. Lett. 122(9), 090601 (2019),  10.1103/PhysRevLett.122.090601,  1810.07170.
  * [13] A. Morvan, T. I. Andersen, X. Mi, C. Neill, A. Petukhov, K. Kechedzhi, D. Abanin, R. Acharya, F. Arute, K. Arya, A. Asfaw, J. Atalaya _et al._ ,  _Formation of robust bound states of interacting microwave photons_ ,  Nature 612, 240 (2022),  10.1038/s41586-022-05348-y,  2206.05254.
  * [14] M. Vanicat, L. Zadnik and T. Prosen,  _Integrable Trotterization: Local Conservation Laws and Boundary Driving_ ,  Phys. Rev. Lett. 121(3), 030606 (2018),  10.1103/PhysRevLett.121.030606,  1712.00431.
  * [15] K. Maruyoshi, T. Okuda, J. W. Pedersen, R. Suzuki, M. Yamazaki and Y. Yoshida,  _Conserved charges in the quantum simulation of integrable spin chains_ ,  J. Phys. A 56(16), 165301 (2023),  10.1088/1751-8121/acc369,  2208.00576.
  * [16] R. J. Baxter,  _Exactly solved models in statistical mechanics_ ,  London: Academic Press Inc (1982).
  * [17] E. Vernier, B. Bertini, G. Giudici and L. Piroli,  _Integrable Digital Quantum Simulation: Generalized Gibbs Ensembles and Trotter Transitions_ ,  Phys. Rev. Lett. 130(26), 260401 (2023),  10.1103/PhysRevLett.130.260401,  2212.06455.
  * [18] I. L. Aleiner,  _Bethe ansatz solutions for certain periodic quantum circuits_ ,  Ann. Phys. 433, 168593 (2021),  10.1016/j.aop.2021.168593,  2107.05715.
  * [19] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean and T. E. O’Brien,  _Quantum error mitigation_ ,  Rev. Mod. Phys. 95(4), 045005 (2023),  10.1103/RevModPhys.95.045005,  2210.00921.
  * [20] W. Hao, R. I. Nepomechie and A. J. Sommese,  _Completeness of solutions of Bethe’s equations_ ,  Phys. Rev E. 88(5), 052113 (2013),  10.1103/PhysRevE.88.052113,  1308.4645.
  * [21] C. Marboe and D. Volin,  _Fast analytic solver of rational Bethe equations_ ,  J. Phys. A 50(20), 204002 (2017),  10.1088/1751-8121/aa6b88,  1608.06504.
  * [22] E. Granet and J. L. Jacobsen,  _On zero-remainder conditions in the Bethe ansatz_ ,  JHEP 03, 178 (2020),  10.1007/JHEP03(2020)178,  1910.07797.
  * [23] Z. Bajnok, E. Granet, J. L. Jacobsen and R. I. Nepomechie,  _On Generalized $Q$-systems_,  JHEP 03, 177 (2020),  10.1007/JHEP03(2020)177,  1910.07805.
  * [24] R. I. Nepomechie,  _Q-systems with boundary parameters_ (2019),  10.1088/1751-8121/ab9386,  1912.12702.
  * [25] R. I. Nepomechie,  _The $A_{m}^{(1)}$ Q-system_,  Mod. Phys. Lett. A 35(31), 2050260 (2020),  10.1142/S0217732320502600,  2003.06823.
  * [26] R. I. Nepomechie,  _Wronskian-type formula for inhomogeneous TQ-equations_ ,  In _International Workshop on Classical and Quantum Integrable Systems_ ,  10.1134/S004057792009007X (2020), 2005.08326.
  * [27] J. Gu, Y. Jiang and M. Sperling,  _Rational $Q$-systems, Higgsing and mirror symmetry_,  SciPost Phys. 14(3), 034 (2023),  10.21468/SciPostPhys.14.3.034,  2208.10047.
  * [28] J. Hou, Y. Jiang and R.-D. Zhu,  _Spin- $s$ Rational $Q$-system_,  SciPost Phys. 16, 113 (2023),  https://scipost.org/10.21468/SciPostPhys.16.4.113,  2303.07640.
  * [29] J. Hou, Y. Jiang and Y. Miao,  _Rational Q-systems at Root of Unity I. Closed Chains_ (2023),  2310.14966.
  * [30] Y. Jiang and Y. Zhang,  _Algebraic geometry and Bethe ansatz. Part I. The quotient ring for BAE_ ,  JHEP 03, 087 (2018),  10.1007/JHEP03(2018)087,  1710.04693.
  * [31] J. Lykke Jacobsen, Y. Jiang and Y. Zhang,  _Torus partition function of the six-vertex model from algebraic geometry_ ,  JHEP 03, 152 (2019),  10.1007/JHEP03(2019)152,  1812.00447.
  * [32] Z. Bajnok, J. L. Jacobsen, Y. Jiang, R. I. Nepomechie and Y. Zhang,  _Cylinder partition function of the 6-vertex model from algebraic geometry_ ,  JHEP 06, 169 (2020),  10.1007/JHEP06(2020)169,  2002.09019.
  * [33] Y. Jiang, R. Wen and Y. Zhang,  _Exact quench dynamics from algebraic geometry_ ,  Phys. Rev. E 108(2), 024128 (2023),  10.1103/PhysRevE.108.024128,  2109.10568.
  * [34] J. Böhm, J. L. Jacobsen, Y. Jiang and Y. Zhang,  _Geometric algebra and algebraic geometry of loop and Potts models_ ,  JHEP 05, 068 (2022),  10.1007/JHEP05(2022)068,  2202.02986.
  * [35] C. Destri and H. De Vega,  _Light-cone lattice approach to fermionic theories in 2D: The massive Thirring model_ ,  Nucl. Phys. B 290, 363 (1987),  10.1016/0550-3213(87)90193-3.
  * [36] C. Destri and H. J. de Vega,  _New thermodynamic bethe ansatz equations without strings_ ,  Phys. Rev. Lett. 69(16), 2313 (1992),  10.1103/PhysRevLett.69.2313.
  * [37] C. Destri and H. J. de Vega,  _Unified approach to Thermodynamic Bethe Ansatz and finite size corrections for lattice models and field theories_ ,  Nucl. Phys. B 438(3), 413 (1995),  DOI: 10.1016/0550-3213(94)00547-R.
  * [38] N. A. Slavnov,  _Algebraic Bethe ansatz_ 10.48550/arXiv.1804.07350,  1804.07350.
  * [39] C. M. Yung and M. T. Batchelor,  _Integrable vertex and loop models on the square lattice with open boundaries via reflection matrices_ ,  Nucl. Phys. B 435, 430 (1995),  10.1016/0550-3213(94)00448-N,  hep-th/9410042.
  * [40] V. E. Korepin,  _Calculation of norms of Bethe wave functions_ ,  Comm. Math. Phys. 86, 391 (1982),  10.1007/BF01212176.
  * [41] A. G. Izergin,  _Partition function of the six-vertex model in a finite volume_ ,  Soviet Physics Doklady 32, 878 (1987).
  * [42] A. G. Izergin, D. A. Coker and V. E. Korepin,  _Determinant formula for the six-vertex model_ ,  J. Phys. A 25(16), 4315 (1992),  10.1088/0305-4470/25/16/010.
  * [43] J. Mossel and J.-S. Caux,  _Relaxation dynamics in the gapped XXZ spin-1/2 chain_ ,  New Journal of Physics 12(5), 055028 (2010),  10.1088/1367-2630/12/5/055028,  1002.3988.
  * [44] J. Salas and A. D. Sokal,  _Transfer matrices and partition function zeros for antiferromagnetic Potts models. 1. General theory and square lattice chromatic polynomial_ ,  J. Statist. Phys. 104, 609 (2001),  10.1023/A:1010376605067,  cond-mat/0004330.
  * [45] K. Fabricius and B. M. McCoy,  _Bethe’s equation is incomplete for the XXZ model at roots of unity_ ,  J. Statist. Phys. 103, 647 (2001),  10.1023/A:1010380116927,  cond-mat/0009279.
  * [46] K. Fabricius and B. M. McCoy,  _Completing Bethe’s equations at roots of unity_ ,  J. Statist. Phys. 104, 573 (2001),  10.1023/A:1010372504158,  cond-mat/0012501.
  * [47] E. Ilievski,  _Popcorn Drude weights from quantum symmetry_ ,  J. Phys. A 55(50), 504005 (2022),  10.1088/1751-8121/acaa77,  2208.01446.
  * [48] D. A. Cox, J. B. Little and D. O’Shea,  _Using Algebraic Geometry_ , vol. 185 of _Graduate Texts in Mathematics_ ,  Springer, first edn.,  10.1007/978-1-4757-6911-1 (1998).
  * [49] W. Decker, G.-M. Greuel, G. Pfister and H. Schönemann,  _Singular 4-3-0 — A computer algebra system for polynomial computations_,  http://www.singular.uni-kl.de (2022).
  * [50] J. Böhm, W. Decker, A. Frühbis-Krüger, F.-J. Pfreundt, M. Rahn and L. Ristau,  _Towards massively parallel computations in algebraic geometry_ (2018), 1808.09727.
  * [51] J. Böhm, H. Xu and Y. Zhang,  _Massively Parallel Methods for Computing Gröbner bases, to appear on arXiv_ (2024).
  * [52] M. Abramowitz, I. A. Stegun and D. Miller,  _Handbook of mathematical functions with formulas, graphs and mathematical tables (national bureau of standards applied mathematics series no. 55)_ ,  Journal of Applied Mechanics 32, 239 (1965).