## Generalized hydrodynamics of integrable quantum circuits

## Friedrich Hübner 1 , Eric Vernier 2 and Lorenzo Piroli 3

- 1 Department of Mathematics, King's College London, Strand, London WC2R 2LS, U.K.
- 2 Laboratoire de Probabilités, Statistique et Modélisation CNRS, Université Paris Cité, Sorbonne Université Paris, France
- 3 Dipartimento di Fisica e Astronomia, Università di Bologna and INFN, Sezione di Bologna, via Irnerio 46, I-40126 Bologna, Italy

## Abstract

Quantum circuits make it possible to simulate the continuous-time dynamics of a manybody Hamiltonian by implementing discrete Trotter steps of duration τ . However, when τ is sufficiently large, the discrete dynamics exhibit qualitative differences compared to the original evolution, potentially displaying novel features and many-body effects. We study an interesting example of this phenomenon, by considering the integrable Trotterization of a prototypical integrable model, the XXZ Heisenberg spin chain. We focus on the well-known bipartition protocol, where two halves of a large system are prepared in different macrostates and suddenly joined together, yielding non-trivial nonequilibrium dynamics. Building upon recent results and adapting the generalized hydrodynamics (GHD) of integrable models, we develop an exact large-scale description of an explicit one-dimensional quantum-circuit setting, where the input left and right qubits are initialized in two distinct product states. We explore the phenomenology predicted by the GHD equations, which depend on the Trotter step and the gate parameters. In some phases of the parameter space, we show that the quantum-circuit large-scale dynamics is qualitatively different compared to the continuous-time evolution. In particular, we find that a single microscopic defect at the junction, such as the addition of a single qubit, may change the nonequilibrium macrostate appearing at late time.

## Contents

| 1   | Introduction             | Introduction                        | 2   |
|-----|--------------------------|-------------------------------------|-----|
| 2   | The model                | The model                           | 4   |
|     | 2.1                      | The integrable quantum circuit      | 4   |
|     | 2.2                      | The Bethe ansatz solution           | 5   |
|     |                          | 2.2.1 The TBA in the gapped regime  | 6   |
|     |                          | 2.2.2 The TBA in the gapless regime | 8   |
| 3   | The bipartition protocol | The bipartition protocol            | 9   |
|     | 3.1                      | The initial states                  | 9   |
|     |                          | 3.1.1                               | 10  |
|     | 3.1.2                    | Gapped phase Gapless phase          | 11  |
|     | 3.2                      | equations                           |     |
|     | GHD                      | GHD                                 | 12  |
| 3.3 | The quasi-momentum and   | the quasi-energy                    | 12  |
|     | 3.3.1                    | Gapped phase                        | 13  |

|            | 3.3.2                                                 | Gapless phase                                         |   14 |
|------------|-------------------------------------------------------|-------------------------------------------------------|------|
| 4          | Results                                               | Results                                               |   14 |
|            | 4.1                                                   | The non-interacting points                            |   15 |
|            | 4.2                                                   | Observable profiles at the interacting points         |   17 |
| 5          | Outlook                                               | Outlook                                               |   18 |
| A          | Numerical solution of the GHD equation                | Numerical solution of the GHD equation                |   19 |
|            | A.1                                                   | Upgraded algorithm                                    |   19 |
| B          | Derivation of the GGE rapidity distribution functions | Derivation of the GGE rapidity distribution functions |   20 |
|            | B.1                                                   | The boundary quantum transfer matrix construction     |   20 |
|            | B.2                                                   | The T-system and Y-system                             |   21 |
|            | B.3                                                   | Truncation of the T- and Y- system at γ = π p + 1     |   22 |
|            | B.4                                                   | The GGE densities at γ = π/ ( p + 1 )                 |   23 |
| C          | Details on the free-fermion points of the model       | Details on the free-fermion points of the model       |   23 |
| References | References                                            | References                                            |   24 |

## 1 Introduction

The Trotter-Suzuki decomposition [ 1,2 ] is at the basis of the idea of quantum simulation [ 3 . ] Up to a controlled error, this decomposition allows us to break down the continuous-time evolution of a many-body system into a series of elementary two-body interactions of duration τ , which could be implemented as quantum gates in a quantum circuit [ 4 . ] The Trotter time τ controls the accuracy and determines the number of gates applied. Unfortunately, external noise limits the number of operations which can be carried out in the noisy-intermediatescale quantum (NISQ) devices [ 5 ] available today. Therefore, since simulating the continuoustime dynamics typically requires application of many quantum gates, full-fledged quantum simulation remains a non-trivial challenge [ 6 . ]

Motivated by this picture, the past few years have witnessed growing interest in genuinely discrete dynamics where the Trotter time τ is not small [ 7,8 . The reason behind this interest ] is two-fold. On the one hand, these models may be easier to implement in current quantum platforms [ 9-13 , since a non-trivial evolution may be observed over a small number of Trotter ] steps. On the other hand, discrete dynamics may exhibit qualitative differences compared to the original Hamiltonian system, with novel phenomena and many-body effects.

For typical Hamiltonians, recent work [ 14 ] has shown that increasing indefinitely the Trotter time τ leads to a sharp transition, see also [ 15-23 : when ] τ is increased beyond a threshold value, approximation errors become uncontrolled at large times. As a consequence, the evolution becomes chaotic, and one expects the rapid onset of local-equilibrium infinite-temperature physics [ 24, 25 . ] This mechanism, however, is not ubiquitous and a different behavior is expected, for instance, in the presence of integrability , where extensively many local conservation laws prevents the onset of quantum chaos [ 26-28 . ]

Integrable quantum circuits providing a Trotterization of integrable Hamiltonians were introduced in Refs. [ 29-31 . ] In the past few years they have received significant attention due to their dynamical features [ 32-35 , their rich underlying mathematical structures ] [ 36-

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

D

S

/

o

L

G

M

3

c

w

7

Y

+

2

P

u

y

&gt;

A

B

H

V

N

8

E

J

U

r

q

9

0

m

k

F

p

Q

z

g

Z

j

I

f

d

W

X

v

5

C

K

T

n

R

O

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

S

v

E

w

9

j

+

G

/

n

N

Z

Q

R

I

&gt;

A

B

H

c

V

8

J

3

U

r

q

L

0

m

k

M

F

2

p

o

z

g

u

D

Y

y

f

d

W

X

7

5

C

K

P

T

O

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

D

S

/

o

L

G

M

3

c

w

7

Y

+

2

P

u

y

&gt;

A

B

H

V

N

8

E

J

U

r

q

9

0

m

k

F

p

Q

z

g

Z

j

I

f

d

W

X

v

5

C

K

T

n

R

O

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

S

v

E

w

9

j

+

G

/

n

N

Z

Q

R

I

&gt;

A

B

H

c

V

8

J

3

U

r

q

L

0

m

k

M

F

2

p

o

z

g

u

D

Y

y

f

d

W

X

7

5

C

K

P

T

O

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

S

v

E

w

9

j

+

G

/

n

N

Z

Q

R

I

&gt;

A

B

H

c

V

8

J

3

U

r

q

L

0

m

k

M

F

2

p

o

z

g

u

D

Y

y

f

d

W

X

7

5

C

K

P

T

O

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

D

S

/

o

L

G

M

3

c

w

7

Y

+

2

P

u

y

&gt;

A

B

H

V

N

8

E

J

U

r

q

9

0

m

k

F

p

Q

z

g

Z

j

I

f

d

W

X

v

5

C

K

T

n

R

O

&lt;latexi sh1\_b64="SvEw9j+G/nNZQRI&gt;ABHcV8J3UrqL0mkMF2pozguDYyfdWX75CKPTO

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

D

S

/

o

L

G

M

3

c

w

7

Y

+

2

P

u

y

&gt;

A

B

H

V

N

8

E

J

U

r

q

9

0

m

k

F

p

Q

z

g

Z

j

I

f

d

W

X

v

5

C

K

T

n

R

O

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-cf30b0b1-6758-4c0c-8de3-a087677003b6.jpg)

-

-

-

-

0 0

.

0 1

.

0 2

.

0 3

.

0 4

.

-

(b)

5

Figure 1: Sketch of the bipartition protocol. The system is initialized by joining at the origin two product states (the figure corresponds to the Néel and anti-Néel states) and evolved by applying a homogeneous quantum circuit, Trotterizing the integrable XXZ continuous-time dynamics ( τ is the duration of the Trotter step). A lightcone spreads ballistically from the junction, inside of which non-trivial dynamics takes place. We study the profiles of local observables arising at late times. The plot shows the example of the staggered spinz magnetization and is obtained by solving the GHD equations, cf. Sec. 4.

39 ] and seminal work demonstrating experimental implementation in current quantum devices [ 40-42 . ] Prominently, Ref. [ 33 ] has shown the existence of Trotter transitions in integrable quantum circuits, manifesting in the correlations of the so-called generalized Gibbs ensemble (GGE) [ 26-28 . ] As GGEs describe the late-time physics of integrable systems, this result is interesting, as it implies that the information on the discretized structure of the dynamics is not washed out at coarse-grained time scales.

In this work, we build upon the results of Ref. [ 33 ] (which focused on homogeneous, translation-invariant settings) and investigate the physics of integrable quantum circuits in more general non-equilibrium situations. We focus in particular on the prototypical model of the XXZ Heisenberg spin chain [ 43 ] and consider the simplest case where the system is prepared in an inhomogeneous state, namely the so-called bipartition protocol [ 44-47 . The ] latter consists in a special type of quantum quench [ 48 , where the two halves of a 1 ] D system are prepared in distinct macrostates, and subsequently left to evolve according to the system unitary dynamics, cf. Fig. 1. While conceptually simple, this problem is hard. For integrable systems, its understanding has required significant theoretical work over the past decade, culminating in the development of the so-called generalized hydrodynamics (GHD) [ 49, 50 , a ] very versatile and powerful theory.

In the bipartition protocol, the GHD predicts that local subsystems at a given distance x from the junction equilibrate to local stationary states, depending only on the ratio x / t , defining a reference frame moving away from the origin. Such stationary states are characterized by space-dependent GGEs, which are quantitatively determined by a solution to the hydrodynamic equations [ 47 . The latter depend on the features of the specific model considered, ] including its particle content and symmetries. For this reason, while GHD provides a very general framework, the qualitative features of the dynamics can vary greatly from model to model, leading to a series of works aimed at exploring the phenomenology of the biparition protocol in different integrable systems [ 49-64 . ]

The goal of this work is to apply the GHD framework to study the bipartition protocol in the XXZ integrable quantum circuit. Based on the known spectral properties of the model, we first derive the GHD equations and study their dependence on the Trotter step and circuit parame-

ζ

0

=

x/t

τ

τ

τ

τ

τ

=0.0

=0.8

=0.4

=1.2

=1.6

5

ters. We consider an explicit protocol where the input left and right qubits are initialized in two distinct product states and evolved by applying a homogeneous quantum circuit. We find that the qualitative features of the discrete dynamics may vary greatly from the continuous-time evolution. Most prominently, we show that a single microscopic defect at the junction, such as the addition of a single qubit, may change the nonequilibrium macrostate appearing at late time. Our work complements Ref. [ 33 , further exploring the interesting dynamical features ] of integrable quantum circuits in many-body settings.

The rest of this work is organized as follows. We begin in Sec. 2 introducing the XXZ quantum circuit and its Bethe ansatz solution. In Sec. 3 we discuss the bipartition protocol and derive the necessary technical ingredients for its analytic description. They include obtaining the GGEs corresponding to different initial states, generalizing the results of Ref. [ 33 , and ] deriving the GHD equations for integrable quantum circuits. Our results are discussed in detail in Sec. 4, where we provide predictions for the system dynamics at different points in the parameter space, by numerically solving the GHD equations. Our conclusions are presented in Sec. 5, while the appendices contain the most technical aspects of our work.

## 2 The model

## 2.1 The integrable quantum circuit

We consider a system of L qubits arranged in a one-dimensional (1 D ) array. We denote by | 0 〉 j , | 1 〉 j the basis states associated with the local Hilbert space at position j . The dynamics consists in the repeated (Floquet) application of the unitary operator

$$U ( \tau ) = U _ { e } ( \tau ) U _ { o } ( \tau ) \,,$$

where

$$U _ { o } ( \tau ) = \prod _ { \substack { n = 1 \\ 1 / 2 } } ^ { L / 2 } V _ { 2 n, 2 n + 1 } ( \tau ),$$

$$U _ { e } ( \tau ) = \prod _ { n = 1 } ^ { \dots } V _ { 2 n - 1, 2 n } ( \tau ) \,,$$

and

$$V _ { n, n + 1 } ( \tau ) = e ^ { - i \frac { \tau } { 4 } \left [ \sigma _ { n } ^ { x } \sigma _ { n + 1 } ^ { x } + \sigma _ { n } ^ { y } \sigma _ { n + 1 } ^ { y } + \Delta ( \sigma _ { n } ^ { x } \sigma _ { n + 1 } ^ { \tau } - 1 ) \right ] } \,,$$

with τ ∈ Ρ , and where periodic boundary conditions are assumed. U ( τ ) defines a brickwork quantum circuit which is a discrete version of the continuous-time evolution driven by the XXZ Heisenberg Hamiltonian [ 29,65 ]

$$H = \frac { 1 } { 4 } \sum _ { j = 1 } ^ { L } \left [ \sigma _ { j } ^ { x } \sigma _ { j + 1 } ^ { x } + \sigma _ { j } ^ { y } \sigma _ { j + 1 } ^ { y } + \Delta ( \sigma _ { j } ^ { z } \sigma _ { j + 1 } ^ { z } - 1 ) \right ] \,,$$

where σ α j are the Pauli matrices acting at position j , with σ α L + 1 = σ α 1 . Indeed, we have

$$e ^ { - i H t } = \lim _ { M \to \infty } U ( t / M ) ^ { M } \,.$$

In other words, for fixed time t , the continuous time evolution is recovered by taking τ → 0 and the circuit depth M →∞ , keeping the product τ M = t constant.

For a given Hamiltonian, the Suzuki-Trotter decomposition is not unique, but the above discretization preserves integrability [ 29, 30 . ] Namely, the Floquet operator (1) features an

extensive number of local and quasilocal conserved operators known as charges [ 29-31 . For ] small τ , they can be thought of as a deformation of those of the Hamiltonian (5): for each charge Qk , with [ Qk , H ] = 0, we have two new operators ˜ Q k ± ( τ ) with [ ˜ Q k ± ( τ ) , U ( τ )] = 0. As pointed out in Refs. [ 29, 31 , the charges ] ˜ Q k ± ( τ ) break the single-site translation symmetry. Instead, the operation leaving the set of conserved charges invariant is the combination of the single-site translation and a change of sign of the staggering parameter x introduced in the next section.

## 2.2 The Bethe ansatz solution

The integrability of the model allows one to obtain the exact spectrum of the Floquet operator U ( τ ) [ 36-38 , which we now review. It is first important to discuss the phase diagram of the ] model. To this end, we introduce the parameters [ 30 ]

$$\gamma = \arccos \left [ \sin ( \Delta \tau / 2 ) / \sin ( \tau / 2 ) \right ] \,,$$

$$x = i \arcsinh \left [ \sin ( \gamma ) \tan ( \tau / 2 ) \right ] \,,$$

and distinguish two cases. First, if γ = i η , with η ∈ Ρ , then x ∈ Ρ . In this case, we say that the model is in the gapped or massive regime. As it will be clear from the subsequent discussion, the name comes from the fact that the spectrum of U ( τ ) can be described in terms of quasiparticles with similar features as those of the XXZ Hamiltonian (5) in the gapped regime | ∆ | &gt; 1. Second, if γ ∈ Ρ then x is purely imaginary. In this case, we say that the model is in the gapless or massless regime. The phase diagram as a function of ∆ and τ , which parametrize the gates (4), can be found in Ref. [ 33 . ] For any ∆ , there is an infinite sequence of transitions as we increase τ . The first phase transition takes place at the value

$$\tau _ { \text{th} } ( \Delta ) = \frac { 2 \pi } { \Delta + 1 } \,.$$

For ∆ &gt; 1 ( ∆ ≤ 1), the system is in a gapped (gapless) phase for τ ≤ τ th ( ∆ ) .

The spectrum of U ( τ ) is organized into sectors labeled by the number M of stable quasiparticles. The latter are magnonic excitations and, accordingly, M is the quantum number associated with the conserved operator ˆ M = P j ( 1 -σ z j ) / 2. In a given M -sector, the eigenstates are parametrized by sets of complex numbers { p j } M j = 1 , satisfying the so-called Bethe equations [ 36 ]

̸

$$\left [ \frac { f _ { x } ^ { + } ( p _ { i } ) } { f _ { x } ^ { - } ( p _ { i } ) } \right ] ^ { \frac { L } { 2 } } = \prod _ { k \neq j } ^ { M } \frac { \sinh \left ( p _ { j } - p _ { k } + i \gamma \right ) } { \sinh \left ( p _ { j } - p _ { k } - i \gamma \right ) } \,,$$

where

$$f _ { x } ^ { \pm } ( p ) = \sinh \left ( p + i \frac { x } { 2 } \pm i \frac { \gamma } { 2 } \right ) \sinh \left ( p - i \frac { x } { 2 } \pm i \frac { \gamma } { 2 } \right ).$$

Each number p j is associated to a stable quasiparticle and is related to the corresponding quasimomentum or rapidity, denoted by λ j . The precise relation between the two depends on the phase of the model: in the gapless regime, we have λ j = p j , while in the gapped regime λ j = i p j . For finite M , it is non-trivial to find all the solutions of the Bethe equations. However, the picture greatly simplifies in the thermodynamic limit L → ∞ , where the spectrum can be described via the so-called Thermodynamics Bethe Ansatz (TBA) formalism [ 66 . ] While initially developed for the study of integrable Hamiltonians, the TBA can be straightforwardly extended to the quantum-circuit setting [ 33 . ]

In essence, the TBA formalism is based on the string hypothesis [ 66 ] which postulates how the solutions to the Bethe equations organize themselves in the complex plane in the limit L →∞ , with the density D = M L / kept fixed. In particular, the string hypothesis states that

the rapidities organize themselves into sets of n elements forming a 'string' of length n , which is interpreted as a bound state of n quasiparticles. Each string features a center, λ , which is the bound-state quasimomentum. Therefore, we can associate to each macrostate a set of functions ρ n ( λ ) : in a large system of size L , L ρ n ( λ ) d λ yields the number of n -quasiparticle bound states with rapidities in the interval [ λ λ , + d λ ] . In order to achieve a complete thermodynamic description, the TBA formalism also introduces the set of hole distribution functions for the n -string centers, denoted by ρ h n ( λ ) . The hole distribution functions are analogous to the distribution of vacancies (i.e., unoccupied states) in the ideal Fermi gas at finite temperature.

The precise form of the rapidity and hole distribution functions depends once again on the phase of the model. We therefore separate the discussion between the two cases in the following. Before proceeding, it is important to mention that the validity of the string hypothesis has been verified in a number of cases, including the present quantum-circuit setting, where the resulting predictions have been tested numerically against independent tensor-network computations [ 33 . ]

## 2.2.1 The TBA in the gapped regime

In the gapped regime, the string length can take any possible integer value, n ≥ 1, while the string centers satisfy λ ∈ [ -π/ 2, π/ 2 . ] By standard techniques [ 66 , one can then take the ] thermodynamic limit of the Bethe equations (9), resulting in an infinite set of coupled integral equations for the distribution functions ρ n ( λ ) , ρ h n ( λ ) . Namely, introducing

$$\rho _ { n } ^ { t } ( \lambda ) = \rho _ { n } ( \lambda ) + \rho _ { n } ^ { h } ( \lambda ) \,,$$

the thermodynamic limit of the Bethe equations read

$$\rho _ { n } ^ { t } ( \lambda ) = a _ { n } ^ { ( x / 2 ) } ( \lambda ) - \sum _ { m = 1 } ^ { \infty } \left ( a _ { n m } * \rho _ { m } \right ) ( \lambda ) \,,$$

where

$$( f * g ) ( \lambda ) \coloneqq \int _ { - \pi / 2 } ^ { \pi / 2 } \mathrm d \mu f ( \mu - \lambda ) g ( \mu ) \,,$$

while we introduced the notation

$$f ^ { ( x ) } ( \lambda ) = ( f ( \lambda + x ) + f ( \lambda - x ) ) / 2 \,,$$

and defined

$$a _ { n m } ( \lambda ) = & ( 1 - \delta _ { n m } ) a _ { | n - m | } ( \lambda ) + 2 a _ { | n - m | + 2 } ( \lambda ) + \dots + 2 a _ { n + m - 2 } ( \lambda ) + a _ { n + m } ( \lambda ) \,, \quad ( 1 5 )$$

with

$$a _ { n } ( \lambda ) = \frac { \sinh ( n \eta ) } { \pi [ \cosh ( n \eta ) - \cos ( 2 \lambda ) ] } \,.$$

Introducing the ratios

$$\eta _ { n } ( \lambda ) = \rho _ { n } ^ { h } ( \lambda ) / \rho _ { n } ( \lambda ) \,,$$

we can also rewrite

$$\rho _ { n } ( \lambda ) [ 1 + \eta _ { n } ( \lambda ) ] = a _ { n } ^ { ( x / 2 ) } ( \lambda ) - \sum _ { m = 1 } ^ { \infty } ( a _ { n m } * \rho _ { m } ) \, ( \lambda ) \,.$$

Note that for x = 0, the above equations recover the standard TBA description of the homogeneous XXZ Heisenberg chain [ 66 . ] Eq. (9) is valid for any eigenstate of the system, and

the dependence in the eigenstate lies in the set of rapidities { p j } . In order to derive Eq. (12) one follows the so-called string hypothesis [ 66 , assuming that in the thermodynamic limit ] the p 's arrange into strings of length n , whose centers are distributed according to ρ n ( λ ) . Eq. (9) then translates into Eq. (12), which relates ρ n ( λ ) to the total density ρ t n ( λ ) . Similar to Eq. (9), Eq. (12) holds for any macrostate, and the dependence in a given macrostate is given by an additional relation between ρ n ( λ ) and ρ t n ( λ ) (in other words, fixing ρ /ρ t n n specifies a macrostate uniquely). While ρ n ( λ ) is sufficient to compute expectation values of conserved quantities, ρ t n ( λ ) is needed to compute other quantities, for instance the entropy of the macrostate ρ n ( λ ) .

As we will see in the next section, the GHD is based on the fact that quasiparticles move ballistically through the system with a given effective velocity v eff n ( λ ) . These velocities are determined the distribution functions ρ n ( λ ) , ρ h n ( λ ) and in particular are obtained as the solution to the system

$$\rho _ { n } ^ { t } ( \lambda ) v _ { n } ^ { \text{eff} } ( \lambda ) & = \frac { 1 } { 2 \pi } E _ { n } ^ { \prime } ( \lambda ) - \sum _ { m = 1 } ^ { \infty } \left ( a _ { n m } * ( \rho _ { m } v _ { m } ^ { \text{eff} } ) \right ) ( \lambda ) \,, & & ( 1 9 ) \\ \dots \, \dots \, \dots \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdots$$

where En ( λ ) will be defined in the following section, cf. Eq. (62). As we will discuss, Eqs. (19) are derived by generalizing the definition given in the Hamiltonian setting in a straightforward way.

In principle, the rapidity distribution functions give us a complete description of the system in the thermodynamic limit. Yet, while simple expressions are known for the expectation value of the local conserved charges, it is generally hard to derive formulas for the correlation functions - a task which, in the XXZ Hamiltonian model, has been carried out only in a few special cases [ 67-70 . In this work, we will restrict ourselves to study three simple observables. ] First, we consider the local averaged magnetization

$$\overline { m } = \frac { \langle \sigma _ { j } ^ { z } \rangle + \langle \sigma _ { j + 1 } ^ { z } \rangle } { 2 } \,.$$

Exploiting the two-site shift-invariance of the model, m coincides with the normalized expectation value of the total spinz magnetization, S z = P j σ z j , which is conserved by the dynamics. In fact, ( 1 -S z ) / 2 counts the number of quasi-particles in each eigenstate, so that we have the simple expression

$$\overline { m } = 1 - 2 \sum _ { n = 1 } ^ { \infty } \int _ { - \pi / 2 } ^ { \pi / 2 } \mathrm d \lambda \, n \rho _ { n } ( \lambda ) \,.$$

The second observable that we consider is the local current associated with the local spinz operator. Its form can be inferred from the continuity equation

$$U ^ { \dagger } ( \tau ) \sigma _ { j } ^ { z } U ( \tau ) - \sigma _ { j } ^ { z } = J _ { j + 1 } - J _ { j } \,,$$

from which we identify

$$J _ { 2 j } = \mathcal { N } _ { x } ^ { - 1 } \left [ 2 ( \sigma _ { 2 j - 1 } ^ { + } \sigma _ { 2 j } ^ { - } - \sigma _ { 2 j - 1 } ^ { - } \sigma _ { 2 j } ^ { + } ) + i \sinh ( i x ) ( \sigma _ { 2 j - 1 } ^ { z } - \sigma _ { 2 j } ^ { z } ) \right ] \,,$$

where

$$\mathcal { N } _ { x } = \frac { i ( 1 + \cosh 2 i x ) } { 2 \sinh i x } \,.$$

Similarly, the expression for J 2 j + 1 is obtained from (23) by shifting all sites by one and changing x to -x , or, equivalently , by simply making the change 2 j -1 → 2 j + 1. The expectation value of

the current associated to a conserved charge is known by recent Bethe Ansatz results [ 71-73 ] and reads

$$\left \langle J _ { 2 j } \right \rangle = - 2 \sum _ { n = 1 } ^ { \infty } \int _ { - \pi / 2 } ^ { \pi / 2 } \mathrm d \lambda n \rho _ { n } ( \lambda ) \nu _ { n } ^ { \text{eff} } ( \lambda ) \,. \quad \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \, \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad. \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \$$

Finally, we consider the one-point function 〈 σ z j 〉 . Since the model breaks single-site translationsymmetry, this expectation value can not be reduced to the value of a conserved charge. In fact, this is a non-trivial one-point function, for which Ref. [ 33 ] derived the formula (expressed in the notation of this paper)

$$\langle \sigma _ { j } ^ { z } \rangle = 1 - 2 \sum _ { n = 1 } ^ { \infty } n \int _ { - \pi / 2 } ^ { \pi / 2 } d \lambda \, \theta _ { n } ( \lambda ) b _ { n } ^ { \text{dr} } ( \lambda ) \,,$$

where we defined b dr n ( λ ) as the solution to the equations

$$b _ { n } ^ { \text{dr} } ( \lambda ) = b _ { n } ( \lambda ) - \sum _ { m } [ a _ { n m } * ( \theta _ { n } ( \lambda ) b _ { m } ^ { \text{dr} } ) ] ( \lambda ) \,,$$

with b n ( λ ) = an ( λ -x / 2 ) and

$$\lambda \, \quad \theta _ { n } ( \lambda ) = \frac { \rho _ { n } ( \lambda ) } { \rho _ { n } ^ { t } ( \lambda ) } \,.$$

Eqs. (21), (25), and (26) will be used to compute the observable profiles arising at late times in the bipartition protocol.

## 2.2.2 The TBA in the gapless regime

The analysis of the gapless regime follows once again the treatment of the Hamiltonian model [ 66 . ] In this case, the string centers satisfy λ ∈ -∞ ∞ ( , ) , while the number and types of allowed strings depend in a non-trivial way on γ . In particular, a simplified structure is achieved at the so-called root-of-unity points, corresponding to γ/π ∈ Θ . In this work, we will restrict to the simplest case where γ takes the form

$$\frac { \gamma } { \pi } = \frac { 1 } { p + 1 } \,,$$

where p ≥ 1 is an integer. In this case, there are Nb = p + 1 different types of strings, which we index by an integer j ∈ { 1, . . . Nb } . Each type is characterized by a number of quasiparticles n j and a parity υ j = ± 1, taking the values

$$\begin{smallmatrix} n _ { j } = j, & v _ { j } = 1, \quad j = 1, 2 \dots, p \,, \\ n _ { p + 1 } = 1, & v _ { p + 1 } = - 1 \,. \end{smallmatrix}$$

Following [ 66 , we also introduce the numbers ]

$$q _ { j } = p + 1 - n _ { j } \quad j = 1, 2 \dots, p \,, & & ( 3 1 ) \\ q _ { p + 1 } = - 1 \,. & & & ( 3 2 )$$

The thermodynamic limit of the Bethe equations (9) then takes the form

$$\text{sign} ( q _ { m } ) \left [ \rho _ { m } ^ { t } ( \lambda ) \right ] = & a _ { m } ^ { ( i x / 2 ) } ( \lambda ) - \sum _ { n = 1 } ^ { N _ { b } } ( a _ { m n } * \rho _ { m } ) \left ( \lambda \right ),$$

where we used the notation (14) and where the convolution is now defined as

$$( f * g ) ( \lambda ) \coloneqq \int _ { - \infty } ^ { \infty } d \mu f ( \mu - \lambda ) g ( \mu ) \,.$$

In addition, we also introduced

$$a _ { j } ( \lambda ) = \frac { v _ { j } } { \pi } \frac { \sin ( \gamma n _ { j } ) } { \cosh ( 2 \lambda ) - v _ { j } \cos ( \gamma n _ { j } ) } \equiv a _ { n _ { j } } ^ { v _ { j } } ( \lambda ) \quad \quad \quad$$

$$a _ { j k } ( \lambda ) = ( 1 - \delta _ { n _ { j } n _ { k } } ) a _ { | n _ { j } - n _ { k } | } ^ { \upsilon _ { j } \upsilon _ { k } } ( \lambda ) + 2 a _ { | n _ { j } - n _ { k } | + 2 } ^ { \upsilon _ { j } \upsilon _ { k } } ( \lambda ) + \dots + 2 a _ { n _ { j } + n _ { k } - 2 } ^ { \upsilon _ { j } \upsilon _ { k } } ( \lambda ) + a _ { n _ { j } + n _ { k } } ^ { \upsilon _ { j } \upsilon _ { k } } ( \lambda ) \,.$$

As in the gapped phase, the quasi-particle rapidity distribution functions allow us to determine the effective velocities, which are obtained as the solution to the equations

$$\text{sign} ( q _ { n } ) \left [ \rho _ { n } ^ { t } ( \lambda ) v _ { n } ^ { \text{eff} } ( \lambda ) \right ] = & \frac { 1 } { 2 \pi } \text{sign} ( q _ { n } ) E _ { n } ^ { \prime } ( \lambda ) - \sum _ { m = 1 } ^ { N _ { b } } \left ( a _ { n m } * ( \rho _ { m } v _ { m } ^ { \text{eff} } ) \right ) ( \lambda ) \,,$$

where E n ′ ( λ ) will be defined later, cf. Eq. (69).

We conclude this section by reporting the expressions, valid in the gapless regime, for the expectation value of the local observables previously introduced. First, the averaged magnetization (20) can be computed via the formula

$$\overline { m } = 1 - 2 \sum _ { j = 1 } ^ { N _ { b } } \int _ { - \infty } ^ { \infty } d \lambda \, n _ { j } \rho _ { j } ( \lambda ) \,,$$

which is analogous to Eq. (21). Second, for the corresponding current we have

$$\langle J _ { 2 n } \rangle = - 2 \sum _ { j = 1 } ^ { N _ { b } } \int _ { - \infty } ^ { \infty } d \lambda n _ { j } \rho _ { j } ( \lambda ) v _ { j } ^ { \text{eff} } ( \lambda ) \,.$$

Finally, the expectation value of the spinz local magnetization reads [ 33 ]

$$\langle \sigma _ { 2 k + m } ^ { \underline { z } } \rangle = 1 - 2 \sum _ { j = 1 } ^ { N _ { s } } \int d \lambda \, n _ { j } \theta _ { j } ( \lambda ) \sigma _ { j } a _ { j } ^ { ( m ) \, \mathrm d r } ( \lambda ) \,,$$

where σ j = sgn ( q j ) , while

$$a _ { j } ^ { ( m ) \, d r } ( \lambda ) = a _ { p } ( \lambda - ( - 1 ) ^ { m \, i x } ) - \sum _ { q = 1 } ^ { N _ { a } } \int d \mu \, a _ { j q } ( \lambda - \mu ) \theta _ { q } ( \mu ) \sigma _ { q } a _ { q } ^ { ( m ) \, d r } ( \mu ),$$

and θ j is defined in Eq. (28).

## 3 The bipartition protocol

## 3.1 The initial states

For systems with traditional Hamiltonian dynamics, most of previous work on bipartition protocols have focused on the setting where two halves of the system are prepared in different equilibrium states. Typically, the latter are chosen as finite-temperature states with different temperatures or chemical potentials [ 44-47, 74-85 . ] In this case, when the two halves are

joined together, the non-trivial dynamics only takes place in a region increasing linearly in time and centered around the junction.

Alternatively, one could consider preparing two halves of the system in non-equilibrium states with a simple structure [ 49 . As discussed in Sec. 1, each half of the system far enough ] from the junction is then expected to relax to a stationary GGE which depends on the initial state. Accordingly, at late time the non-trivial dynamics again takes place only in a region increasing linearly in time and centered around the junction.

In the quantum-circuit setting under consideration, the Hamiltonian (5) is not conserved by the discrete time evolution, and the corresponding Gibbs states are not equilibrium states. Therefore, the most natural choice is to initialize the two halves of the system into different nonequilibrium states, which we choose as product states. This setting is natural in light of benchmark implementation in quantum devices, since product states are the easiest to be prepared in digital platforms.

Our analysis relies on the possibility to characterize the GGE associated to the chosen initial states. As shown in Ref. [ 33 , this problem can be solved for a large family of states ] called integrable, generalizing the constructions of the Hamiltonian case [ 86-88 . The family ] of integrable states includes all two-site shift-invariant product states, allowing us to study a wide class of bipartition protocols. For concreteness, we will focus on the following states:

- • the Néel state

$$| \Psi _ { N } \rangle = | 0 \rangle \otimes | 1 \rangle \otimes \cdots | 0 \rangle \otimes | 1 \rangle \,,$$

- • the anti-Néel state,

$$| \Psi _ { \text{AN} } \rangle = | 1 \rangle \otimes | 0 \rangle \otimes \cdots | 1 \rangle \otimes | 0 \rangle \,,$$

- • the Majumdar-Ghosh or dimer state,

$$| \Psi _ { \text{D} } \rangle = \left ( \frac { | 0 1 \rangle - | 1 0 \rangle } { \sqrt { 2 } } \right ) ^ { \otimes L / 2 } \,.$$

In the following, we provide analytic formulas for the rapidity distribution functions of GGEs corresponding to these states. As usual, we give a separate discussion for the gapped and gapless regime. As the derivation of these results is quite technical, we report it in Appendix B, to which the reader is referred to for further details.

## 3.1.1 Gapped phase

In the gapped phase, the GGE is characterized by an infinite set of functions { η j } j ≥ 1 , defined in Eq. (17). For the states under consideration, all those functions follow from the knowledge of η 0 = 0 and η 1 through iterated applications of the following 'Y-system' relation

$$\eta _ { j } \left ( \lambda + i \frac { \gamma } { 2 } \right ) \eta _ { j } \left ( \lambda - i \frac { \gamma } { 2 } \right ) = & \left [ 1 + \eta _ { j + 1 } \left ( \lambda \right ) \right ] \left [ 1 + \eta _ { j - 1 } \left ( \lambda \right ) \right ],$$

where we recall that γ = i η is purely imaginary throughout the gapped phase.

The function η 1 , in turn, is given by

$$1 + \eta _ { 1 } ( \lambda ) = \left ( 1 + \mathfrak { a } ( i \lambda - \gamma / 2 ) \right ) \left ( 1 + 1 / \mathfrak { a } ( i \lambda + \gamma / 2 ) \right ) \,.$$

For the Néel and anti-Néel state, we find

$$\mathfrak { a } ( u ) = \frac { \sin ( 2 u + \gamma ) } { \sin ( 2 u - \gamma ) } \frac { \sin \left ( u \pm \frac { x } { 2 } \right ) } { \sin \left ( u \mp \frac { x } { 2 } \right ) } \frac { \sin \left ( u - \gamma \pm \frac { x } { 2 } \right ) } { \sin \left ( u + \gamma \mp \frac { x } { 2 } \right ) }$$

Figure 2: GGE rapidity distribution functions ρ n ( λ ) for the Néel initial state. (a): ρ 1 ( λ ) in the gapped phase ∆ = 2, and different values of τ . (b), (c): same plot for ρ 1 ( λ ) and ρ 3 ( λ ) , respectively, in the gapless regimes and different values of τ . The value of ∆ is chosen so that γ = π/ 3, cf. Eq. (7a).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-808cef96-cb07-4727-ba67-76863d97d08c.jpg)

while for the dimer state

$$\mathfrak { a } ( u ) = \frac { \tan ( u + \gamma ) } { \tan ( u - \gamma ) } \,, \quad \ \quad \ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \,. \ \. \ \. \ \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \$$

see Appendix B. Note that the fact that the GGEs arising from the Néel and anti-Néel states display the same rapidity distribution functions in the gapped phase implies that they can not be distinguished by any local correlation function.

Given the functions η n ( λ ) , the rapidity distribution functions are obtained solving Eqs. (18), which can be done by truncating the infinite system and using standard iterative methods. An example of a numerical solution is displayed in Fig. 2. In the limit x = 0, the rapidity distribution functions above recover the well-known results for the Hamiltonian case [ 68,89-94 . ] Note that, while the dependence in the parameter x disappears in the dimer case, the GGE densities ρ n depend on it as x appears in Eq. (18).

## 3.1.2 Gapless phase

As stated previously, the GGE description in the gapless phase depends in a non-trivial way on the specific value of γ . Following Section 2.2, we shall limit ourselves to the simplest case where γ takes the form (29). The GGE is then characterized by a family of functions { η j } j = 1,... p + 1 , which are obtained from a set of auxiliary functions { Yj } j = 1,... p , K , b whose construction is detailed in Appendix B. For concreteness, we give below the explicit rapidity distributions for p = 2, 3. For p = 2 and the Néel state, we have

$$\eta _ { 1 } ( \lambda ) = \eta _ { 2 } ( \lambda ) = \frac { \cosh ( 2 \lambda + i x ) } { \cosh ( 2 \lambda - i x ) }$$

while for the dimer,

$$\eta _ { 1 } ( \lambda ) = 1 / \eta _ { 2 } ( \lambda ) = ( \tanh \lambda ) ^ { 2 } \,.$$

For p = 3, we have, for the Néel state,

$$- \, \wp, \, \intertext { a v e, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i ou, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o u, i o. \, \intertext { a v e }$$

and for the dimer,

$$\eta _ { 1 } ( \lambda ) & = 3 \tanh ^ { 2 } ( \lambda ) \,, \\ \eta _ { 2 } ( \lambda ) & = 2 - \frac { 6 } { 1 + 2 \cosh ( 2 \lambda ) } \,, \\ \eta _ { 3 } ( \lambda ) & = \frac { 1 } { 8 } \left ( 4 + \frac { 3 } { \sinh ( \lambda ) ^ { 2 } } \right ) \,. \\ \intertext { \ } \text{ be plugged into the Bethe equations (33) to obtain the corresponding}$$

These expressions can be plugged into the Bethe equations (33) to obtain the corresponding rapidity distribution functions, see Fig. 2 for an example.

## 3.2 GHD equations

The theory of GHD describes the system in the Euler scaling limit (large times and large system sizes). In the Hamiltonian setting, its formulation is based on the TBA description of the model, from which one can read off three necessary ingredients: The quasi-momentum Pn ( λ ) and the energy En ( λ ) of each quasi-particle type, and the two-particle scattering phase shifts ϕ nm ( λ µ , ) between two strings of lengths n , m and centers λ µ , , respectively.

In our work, we face the problem of generalizing the standard GHD framework to the quantum-circuit setting. To this end, it is useful to view the two-step evolution operator as generated by the Floquet Hamiltonian H F , defined via

$$U ( \tau ) = e ^ { - i \tau H _ { F } } \,.$$

Compared to the continuous-time setting, the Floquet Hamiltonian is not local. However, the corresponding charges are local [ 29 , and we can repeat the reasoning leading to the derivation ] of the GHD equations in the continuous-time setting [ 49,50 . Note also that all eigenvectors of ] U ( τ ) (which are given by the Bethe ansatz) are also eigenvectors of H F and the corresponding eigenvalues En ( λ ) (the quasi-energies) are unique up to integer multiples of 2 π . However, this additive constant is immaterial, since the GHD equation only depends on the derivatives E n ′ ( λ ) .

In summary, we can apply the standard GHD prescriptions to the quantum circuit by replacing the energy of a quasi-particle by its quasi-energy En ( λ ) . Using these definitions, the GHD equations take the standard form [ 49,50 ]

$$\partial _ { t } \rho _ { n } ( t, x, \lambda ) + \partial _ { x } ( v _ { n } ^ { \text{eff} } ( t, x, \lambda ) \rho _ { n } ( t, x, \lambda ) ) = 0,$$

where the effective velocity of string with length n and rapidity λ satisfies the following selfconsistency equation depending on ρ n ( t , x , λ )

$$v _ { n } ^ { \text{eff} } ( t, x, \lambda ) = \frac { E ^ { \prime } _ { n } ( \lambda ) } { P ^ { \prime } _ { n } ( \lambda ) } + \sum _ { m } \int d \mu \, \frac { \varphi _ { m m } ( \lambda, \mu ) } { P ^ { \prime } _ { n } ( \lambda ) } \rho _ { m } ( t, x, \mu ) ( v _ { m } ^ { \text{eff} } ( t, x, \mu ) - v _ { n } ^ { \text{eff} } ( t, x, \mu ) ) \,,$$

which can be rewritten as

$$\rho _ { n } ^ { t } ( \lambda ) v _ { n } ^ { \text{eff} } ( \lambda ) = \frac { 1 } { 2 \pi } E _ { n } ^ { \prime } ( \lambda ) + \sum _ { m } \int d \mu \, \varphi _ { n m } ( \lambda, \mu ) \rho _ { m } ( \mu ) v _ { m } ^ { \text{eff} } ( \mu ).$$

As we will discuss in the next subsection, the form of the scattering phase depends on the phase of the model and it can be read off by comparing (56) with Eqs. (19) and (37), respectively.

## 3.3 The quasi-momentum and the quasi-energy

We finally detail the identification of the quasi-momenta and quasi-energies, separating the discussion for the gapped and gapless phases.

## 3.3.1 Gapped phase

As discussed in Sec. 2, the eigenvalues of U ( τ ) depend on the solutions { p j } to the Bethe equation (9). Denoting by u ( { p j } ) the eigenvalue of U ( τ ) corresponding to a set { p j } M j = 1 , we have

$$u ( \{ p _ { j } \} ) \ = \ \prod _ { j = 1 } ^ { M } g ^ { + } ( p _ { j } ) g ^ { - } ( p _ { j } ) \,,$$

where

$$g ^ { + } ( p ) = \frac { \sinh \left ( p + i \frac { x } { 2 } - i \frac { \gamma } { 2 } \right ) } { \sinh \left ( p + i \frac { x } { 2 } + i \frac { \gamma } { 2 } \right ) } \,, \\ \dots \, \tau \quad. \, \tilde { \ } \quad. \, \tilde { \ } \,. \, \nu \, \, \mathbb { \ } }$$

$$g ^ { - } ( p ) = \frac { \overset { \vee } { \sinh \left ( p - i \frac { x } { 2 } + i \frac { \gamma } { 2 } \right ) } } { \sinh \left ( p - i \frac { x } { 2 } - i \frac { \gamma } { 2 } \right ) } \,.$$

This result follows from the exact solution of the model, see e.g. Ref. [ 33,36 . From (57), the ] quasienergies are expressed as a sum over the individual rapidities, namely

$$E ( \{ p _ { j } \} ) = \sum _ { j } \epsilon ( p _ { j } ) \,,$$

where

$$\epsilon ( p ) = - \frac { i } { \tau } \log \frac { \sinh ( p + i \frac { x } { 2 } - i \frac { \gamma } { 2 } ) } { \sinh ( p + i \frac { x } { 2 } + i \frac { \gamma } { 2 } ) } \frac { \sinh ( p - i \frac { x } { 2 } + i \frac { \gamma } { 2 } ) } { \sinh ( p - i \frac { x } { 2 } - i \frac { \gamma } { 2 } ) } \,.$$

Note that we have divided by τ , consistent with the definition of H F in Eq. (53). We can now write the quasienergy of a n -string by summing over the n rapidities composing it. Recalling that in the gapped phase the physical rapidities are λ j = i p j , we have

$$E _ { n } ( \lambda ) = \sum _ { j = - \frac { n - 1 } { 2 } } ^ { \frac { n - 1 } { 2 } } \varepsilon ( i ( \lambda + j \eta ) ) = - \frac { i } { \tau } \log \left [ \frac { \sin \left ( \lambda - \frac { x } { 2 } + i n \frac { \eta } { 2 } \right ) \sin \left ( \lambda + \frac { x } { 2 } - i n \frac { \eta } { 2 } \right ) } { \sin \left ( \lambda - \frac { x } { 2 } - i n \frac { \eta } { 2 } \right ) \sin \left ( \lambda + \frac { x } { 2 } + i n \frac { \eta } { 2 } \right ) } \right ].$$

Finally, taking the derivative with respect to λ , we arrive at the result

$$E _ { n } ^ { \prime } ( \lambda ) = \frac { - \sinh ( n \eta ) / \tau } { \cosh ^ { 2 } ( \frac { n \eta } { 2 } ) \sin ^ { 2 } ( \lambda - \frac { x } { 2 } ) + \sinh ^ { 2 } ( \frac { n \eta } { 2 } ) \cos ^ { 2 } ( \lambda - \frac { x } { 2 } ) } - ( x \to - x ) \,.$$

Note that, in the limit τ → 0, the above formula recovers the known expression for the continuous time evolution under the XXZ Hamiltonian

$$E _ { n } ^ { \prime } ( \lambda ) \rightarrow - \pi \sinh ( \eta ) a _ { n } ^ { \prime } ( \lambda ) \,.$$

The quasi-momenta can be determined from the Bethe equations (18). In summary, in the gapped phase we have

$$P _ { n } ^ { \prime } ( \lambda ) = 2 \pi a _ { n } ^ { ( x / 2 ) } \,,$$

$$E _ { n } ^ { \prime } ( \lambda ) = 2 \pi ( a _ { n } ( \lambda + x / 2 ) - a _ { n } ( \lambda - x / 2 ) ) / \tau \,,$$

while the scattering phase of two strings reads

$$\varphi _ { n m } ( \lambda, \mu ) = - a _ { n m } ( \lambda - \mu ) \,.$$

Figure 3: Late time profiles after bipartition protocols at the non-interacting points, τ = π/ 2, ∆ = 4.0. (a): the left and right states are the Néel and dimer states, respectively . The plot shows the profiles for the local magnetization, mj =( 〈 σ z j 〉 + 〈 σ z j + 1 〉 ) / 2, and the corresponding current 〈 J j 〉 . Gray lines are the solution to the GHD equations, while blue and orange lines correspond to the exact simulation of the dynamics using the Gaussian formalism. The latter displays small ripples, which are finite-time effects. (b): same plot, for the protocol in which the left and right states are the Néel and anti-Néel states. All profiles of the Gaussian numerical computations are obtained at time t = 256 τ and were carried out on a sufficiently large system of size L = 4096 (note that, due to the strict lightcone on a quantum circuit, this is equivalent to simulating the infinite system).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-c8ff4219-174e-4c85-a53c-b7c779231663.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-3578110c-ceb6-4f1c-b17e-9d31192db2ce.jpg)

## 3.3.2 Gapless phase

The expressions in the gapless phase can be obtained from those in the gapped regime by replacing λ →-i λ and η → i γ . Taking into account the sign sgn ( q n ) appearing in the Bethe equations, we recover the standard formulation of GHD by identifying

$$P _ { n } ^ { \prime } ( \lambda ) = 2 \pi \, \text{sgn} ( q _ { n } ) a _ { n } ^ { ( i x / 2 ) } \,,$$

$$E _ { n } ^ { \prime } ( \lambda ) = 2 \pi \frac { \text{sgn} ( q _ { n } ) } { \tau } \left [ a _ { n } \left ( \lambda - i \frac { x } { 2 } \right ) - a _ { n } \left ( \lambda + i \frac { x } { 2 } \right ) \right ] \,,$$

while the scattering phase is

$$\varphi _ { n m } ( \lambda, \mu ) = - \text{sgn} ( q _ { n } ) a _ { n m } ( \lambda - \mu ) \,.$$

## 4 Results

Having derived the GHD equations for our model, we now proceed to study them. First, following the standard theory [ 49, 50 , we recall that the GHD equations can be written as ] differential equations in terms of the rescaled variable

$$\zeta = x / t \,,$$

namely

$$\zeta \partial _ { \zeta } \rho _ { n } ( \zeta, \lambda ) = \partial _ { \zeta } ( v _ { n } ^ { \text{eff} } ( \zeta, \lambda ) \rho _ { n } ( \zeta, \lambda ) ).$$

Eq. (72) describes how the local, space-time dependent GGE changes as we vary the velocity ζ of the reference frame. In the bipartition protocols, this equation is supplemented by the

boundary conditions

$$\lim _ { \zeta \to \infty } \rho _ { n } ( \zeta, \lambda ) = \rho _ { n } ^ { R } ( \lambda ) \,,$$

$$\lim _ { \zeta \to - \infty } \rho _ { n } ( \zeta, \lambda ) = \rho _ { n } ^ { L } ( \lambda ) \,,$$

where ρ L n ( λ ) , ρ R n ( λ ) are the GGE distribution functions associated with the left and right initial states. For the product states considered in this work, these functions were computed in Sec. 3. The above boundary conditions encode the fact that, away from the junction point, the system will locally equilibrate to a homogeneous GGE, determined by the initial states.

Solving (72) with the boundary conditions (73) can be done with standard iterative methods, see e.g. Ref. [ 95 . ] In this work, we use a different numerical approach which was recently introduced in Ref. [ 96 ] and briefly explained in Appendix A. In the rest of this section we present the final results of our numerical computations. Note that the solution to (72) gives us a family of rapidity distribution functions { ρ n ( ζ λ , ) } ζ . We use the latter to compute the space-time profiles of local observables, using the formulas presented in Sec. 3.

## 4.1 The non-interacting points

We begin by discussing the non-interacting points of the model, where we can test our predictions against independent, numerically exact results. These points correspond to either one of the following choices of the circuit parameters [ 33 : ]

̸

- (i) ∆ = 0 and τ = 2 π / n ∆ , with n ∈ Ζ ;
- (ii) ∆ = 0 and τ ∈ Ρ + .

At these special values of the parameter space, we can make use of the Jordan-Wigner transformation

$$c _ { n } ^ { \dagger } = \frac { \sigma _ { n } ^ { x } + i \sigma _ { n } ^ { y } } { 2 } \prod _ { j = 1 } ^ { n - 1 } \sigma _ { j } ^ { z }, \quad c _ { n } = \frac { \sigma _ { n } ^ { x } - i \sigma _ { n } ^ { y } } { 2 } \prod _ { j = 1 } ^ { n - 1 } \sigma _ { j } ^ { z } \,,$$

to rewrite

$$U _ { o } ( \tau ) = \exp \left [ - i \frac { \tau } { 2 } \sum _ { n = 1 } ^ { L / 2 } \left ( c _ { 2 n } ^ { \dagger } c _ { 2 n + 1 } + c _ { 2 n + 1 } ^ { \dagger } c _ { 2 n } \right ) \right ], \quad \quad \quad$$

$$U _ { e } ( \tau ) = \exp \left [ - i \frac { \tau } { 2 } \sum _ { n = 1 } ^ { L / 2 } \left ( c _ { 2 n + 1 } ^ { \dagger } c _ { 2 n + 2 } + c _ { 2 n + 2 } ^ { \dagger } c _ { 2 n + 1 } \right ) \right ],$$

where { c n } , { c † n } satisfy the canonical anti-commutation relations.

Eqs. (75) and (76) make it explicit that the operators Ue ( τ ) , Uo ( τ ) are written as the exponential of Hamiltonians which are quadratic in the fermion operators and thus noninteracting. As we explain in Appendix C, the corresponding dynamics can then be computed in a numerically-exact and efficient way by means of the so-called Gaussian formalism [ 97 , ] allowing us to produce data for large system sizes and simulation times. These computations serve as an important benchmark for the GHD equations and our numerical solution to them. As we discuss in the following, the GHD predictions are always found to be in perfect agreement with the exact numerics, up to controlled and expected finite-size and finite-time effects.

In Fig. 3(a), we first show data for the protocol where the left and right halves of the system are initialized in the Néel and dimer state respectively. The plot shows the profiles

Figure 4: Late time profiles after the bipartition protocol in which the left and right states are the Néel and dimer states, respectively. The plots correspond to (a) the local staggered magnetization ( 〈 σ z j 〉-〈 σ z j + 1 〉 ) / 2 and (b) the averaged local magnetization mj = ( 〈 σ z j 〉 + 〈 σ z j + 1 〉 ) / 2. We set τ = 1.2, while the value of ∆ is chosen in such a way that γ = π/ 3.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-d7a3e5be-bfa6-4ab8-9e72-cada2bdb8f42.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-ca736a77-a547-4439-a0bb-6c70d33aad43.jpg)

of the local magnetization and the corresponding current. We see that the profiles are flat outside the 'lightcone' [ -ζ max , ζ max ] , which is an expected feature: indeed ζ max corresponds to the maximum velocity at which quasi-particles propagate. Taking reference frames moving at a velocity with larger absolute value, the system is described by the GGEs associated by the initial left / right states.

The magnetization is vanishing at ζ = ±∞ , but the profiles is non-trivial inside of the lightcone, signaling the emergence of non-trivial, non-equilibrium steady states. We also note that the magnetization and current profiles display no points of non-analiticities inside of the lightcone. This is expected, since at the non-interacting points we have a single quasi-particle string type [ 47 . ]

The features described above are qualitatively similar to what one could have expected based on the GHD in the XXZ Heisenberg model [ 49,53 . Yet, away from the limit of small Trot-] ter step, one may expect the emergence of qualitative differences compared to the continuoustime evolution [ 33 . We can see that this is indeed the case even at the non-interacting point ] of the model with non-zero τ . To this end, we consider the bipartition protocol where the left and right states are prepared in a Néel and anti-Néel state, respectively. Note that the resulting initial state can be thought of as a homogeneous Néel state in which a defect, in the form of a | 0 〉 qubit, is added in the center of the chain.

This is a peculiar situation because the Néel and anti-Néel states are related by a translation of a single site. In models in which all charges satisfy single-site translation symmetry, such as the XXZ Heisenberg chain, both states equilibrate to the same GGE [ 98 . Accordingly, after ] an initial transient time, there would be no non-trivial dynamics at large space-time scales. This is intuitive when viewing the initial state as created by the addition of a single defect in an homogeneous configuration. Indeed, a local initial perturbation does not typically lead to a change in the macroscopic state of the system, although exceptions exist, see e.g. [ ? , ? ? , , 99-101 . ] On the other hand, the Trotterized XXZ circuit is invariant under translation by two sites, which allows for the possibility that the Néel and anti-Néel states give rise to different GGEs. It is worth mentioning that one could construct local integrable Hamiltonians breaking the single-site translation symmetry (but preserving two-site shift invariance) for which a similar behavior is observed. However, while this construction would appear rather ad hoc in the continuous-time limit, here it appears naturally as a consequence of increasing the Trotter step. The ability to change the global late-time behavior by a single defect in a circuit setting is also of obvious interest for benchmark experiments in current digital devices.

Figure 5: Late time profiles after the bipartition protocol in which the left and right states are the Néel and anti-Néel states. The plots show (a) the local staggered magnetization ( 〈 σ z j 〉 - 〈 σ z j + 1 〉 ) / 2, (b) the averaged local magnetization mj =( 〈 σ z j 〉 + 〈 σ z j + 1 〉 ) / 2, and (c) the current 〈 J j 〉 . Each curve corresponds to a different τ , while the value of ∆ is always chosen in such a way that γ = π/ 3.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00474v5-figure-d9d33619-ba6b-41a5-bd07-be02fbf1fdfb.jpg)

Indeed, approximating the continuous-time evolution of a two-site shift-invariant Hamiltonian would require implementing very deep circuits, even to reach short times, which is unfeasible on current devices. Instead, existing experiments [ 40-42 ] have already shown that integrable circuits can be implemented keeping sufficient coherence over a number of discrete steps that is enough to appreciate the emergence of late-time features.

Fig. 3(b) shows our quantitative results for this bipartition protocol in the non-interacting case. The profiles of the local magnetization and current demonstrate the emergence of nontrivial dynamics even at macroscopic scales. Interestingly, we see that there is an accumulation of the magnetization charge inside of the lightcone, with the current changing sign at the left and right ends. This does not violate the global conservation of the z -magnetization in the system, as the excess charge propagates to infinity. It is natural to ask whether this peculiar behaviour is a special property of the non-interacting points. We see in the next section that this is not the case, and that non-trivial dynamics in this bipartition protocol survives the presence of interactions.

## 4.2 Observable profiles at the interacting points

We finally present our results for different bipartition protocols at generic points of the parameter space.

First, it is important to state that, for the initial states considered, the magnetization and current profiles are trivial in the gapped regime. This is because both the Néel and dimer states have vanishing magnetization, and in the gapped regime there are no additional charges breaking the spin-flip symmetry [ 33 . We have numerically verified this prediction by solving ] the GHD equations and computing the magnetization and current profiles. In all cases, we have found that the profiles of these observables are flat, up to the chosen numerical accuracy. Note, that the solutions { ρ n ( ζ λ , ) } to the GHD equations for the Néel-dimer bipartition protocol do have a non-trivial dependence on ζ , even in the gapped regime. Conversely, the GGEs associated with the Néel and anti-Néel states coincide, and there is no non-trivial dynamics taking place in the corresponding bipartition protocol.

Next, we move on to discuss the gapless phase. In Fig. 4, we study the protocol where the left and right states are initialized in the dimer and Néel states, respectively . We plot the profiles of the staggered and averaged magnetization, respectively, both showing non-trivial dynamics. Compared to the non-interacting circuit, we see additional points inside of the lightcone where

the profiles are not smooth. Similarly to the GHD profiles in the XXZ Heisenberg chain [ 49,53 , ] these points are associated to the maximum velocities of distinct quasi-particle types (namely, of strings with increasing length n ). We also see that both observables display non-monotonic and highly irregular profiles. While the figure shows data for one value of τ and ∆ , we have verified that qualitatively similar features arise for different choices of the circuit parameters.

Finally, we discuss the bipartition protocol where the two halves of the system are prepared in the Néel and anti-Néel state, respectively. We report an example of our results in Fig. 5, where we plot the profiles of the staggered and averaged magnetization, and the current (again in the gapless phase). As anticipated, all plots show non-trivial profiles and, once again, we see non-analyticity points inside of the lightcone, although they are less evident than in the Néeldimer bipartition protocol. As in the non-interacting case, we see that the dynamics causes an accumulation of magnetization at the origin [ Fig. 5(b) , with the current taking opposite ] values at the left and right ends of the lightcone.

## 5 Outlook

In this work, we have developed and studied the GHD of integrable quantum circuits. Our main finding is that, for generic choices of the circuit parameters, the large-scale dynamics displays qualitative differences compared to the continuous-time evolution of the corresponding XXZ Hamiltonian. In particular, considering quenches from product states, we predict peculiar features which can be detected by looking at the late-time dynamics of simple observables, such as the local magnetization.

Our work opens several directions for future research. First, it would be important to provide predictions for more general observables and correlations than those considered in this work. Technically, this problem requires the derivation of new formulas for the expectation value of local observables in stationary states described by a given rapidity distribution function. While this is a hard problem, we expect that one could generalize the results available for the Hamiltonian model [ 67-70 , extending the computation presented in Ref. ] [ 33 ] for the simplest case of one-point functions.

Along similar lines, it would be interesting to also consider the entanglement dynamics in integrable quantum circuits. In the Hamiltonian setting, building upon the results for the evolution of bipartite entanglement in homogeneous systems [ 102 , recent work has shown ] how the GHD formalism is naturally suited to describe the entanglement dynamics in inhomogeneous setting [ 47,55 . We expect that the results of Refs. ] [ 102 ] apply directly to integrable quantum circuits, potentially allowing us to investigate the role of discrete space-time in the large-scale dynamics of bipartite entanglement.

Finally, it would be important to analyze the feasibility of observing some of the features predicted by our work in benchmark experiments in existing quantum platforms [ 40-42 . To ] this end, it would be necessary to study how finite-size corrections affect our predictions. In addition, one would also need to analyze the effect of noise, possibly taking into account recent strategies for quantum error-mitigation [ 103-112 . ] We leave these interesting questions for future work.

## Acknowledgements

We acknowledge helpful discussions with Bruno Bertini and especially Benjamin Doyon for early collaboration on this project.

Funding information FH acknowledges funding from the faculty of Natural, Mathematical &amp; Engineering Sciences at King's College London. The CREATE cluster [ 113 ] was used for numerical simulations. This work was co-funded by the European Union (ERC, QUANTHEM, 101114881). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## A Numerical solution of the GHD equation

We solve the GHD equation numerically using a recently introduced fixed point equation [ 96 . ] Given the states on the left ( -) and the right side ( + ) either in terms of their quasi-particle densities ρ ± ; i ( λ ) or equivalently as occupation functions θ ± ; i ( λ ) , we can compute

$$P ^ { \prime \text{dr} } _ { \pm ; i } ( \lambda ) = P ^ { \prime } _ { i } ( \lambda ) + \sum _ { j } \int \text{d} \mu \, \varphi _ { i j } ( \lambda - \mu ) \rho _ { \pm ; j } ( \mu ).$$

Then define:

$$K _ { 0 ; i } ( x, \lambda ) = P _ { \ s g n ( x ), i } ^ { \prime \, \text{dr} } ( \lambda ) x$$

$$\hat { \theta } _ { 0 ; i } ( k, p ) = \theta _ { \text{sgn} ( k ) ; i } ( \lambda )$$

$$\hat { N } _ { 0 ; i } ( k, p ) = \theta _ { \text{sgn} ( k ) ; i } ( \lambda ) k$$

As derived in Ref. [ 96 , ] x → K 0; i ( x , λ ) can be interpreted as a change of spatial coordinates system that trivializes the GHD equation, i.e. ∂ θ t ˆ i ( k , λ ) = 0. To obtain the solution at spacetime point ( t , x ) , we thus only need to find Ki ( t , x , λ ) , which satisfies the following integralfixed-point equation [ 96 : ]

$$& K _ { i } ( t, x, \lambda ) = P ^ { \prime } _ { i } ( \lambda ) x - E ^ { \prime } _ { i } ( \lambda ) t + \hat { \mathbf T } \hat { N } _ { 0 ; j } ( K _ { j } ( t, x, \mu ), \mu ) \\ & = P ^ { \prime } _ { i } ( \lambda ) x - E ^ { \prime } _ { i } ( \lambda ) t + \hat { \mathbf T } [ \theta _ { \text{sgn} ( K _ { j } ( t, x, \mu ) ) ; j } ( \mu ) K _ { j } ( t, x, \mu ) ].$$

Here ˆ T f i ( λ ) = P R j d µϕ i j ( λ µ , ) f j ( µ ) . Once the solution Ki ( t , x , λ ) is obtained we can compute the solution to the GHD equation as follows:

$$\theta _ { i } ( t, x, \lambda ) = \hat { \theta } _ { i } ( K _ { i } ( t, x, \lambda ), \lambda ) = \theta _ { \text{sgn} ( K _ { i } ( t, x, \lambda ) ) ; i } ( \lambda ).$$

We solve Eq. (A.5) by iteration with initial guess P i ′ ( λ ) x -E i ′ ( λ ) t , which converges in the gapless phase, where only finitely many strings are present.

## A.1 Upgraded algorithm

As it turns out for the Néel and Dimer state in the gapped phase this iterative scheme does not converge. The same holds for the dressing equation (which can be viewed as a linearized version of (A.5)). The iteration scheme does not converge because ni ( λ ) does not decay sufficiently fast for large strings i → ∞ : Consider the first term of the iteration procedure to compute f dr i ( λ )

$$f _ { i } ^ { \text{dr} } ( \lambda ) = f _ { i } ( \lambda ) + \hat { \mathfrak T } \theta _ { j } ( \mu ) f _ { j } ( \mu ) + \dots$$

Since ni ( λ ) does not decay fast enough, already the first term ˆ T n j ( µ ) f j ( µ ) diverges, unless f i ( λ ) decays sufficiently fast (which it does not for instance for f i ( λ ) = P ′ i ( λ ) ). We can however slightly modify the iteration procedure to obtain a convergent fixed-point iteration.

Let us first explain the idea starting with the dressing equation:

$$f _ { i } ^ { \text{dr} } ( \lambda ) = f _ { i } ( \lambda ) + \hat { \mathbf T } \theta _ { j } ( \mu ) f _ { j } ^ { \text{dr} } ( \mu ).$$

We can rewrite this as follows: Choose any functions g i ( λ ) and mi ( λ ) and define

$$b _ { i } ( \lambda ) = g _ { i } ( \lambda ) - \hat { \mathfrak T } m _ { j } ( \mu ) g _ { j } ( \mu ).$$

Now write f dr i ( λ ) = wi ( λ ) g i ( λ ) and observe:

$$w _ { i } ( \lambda ) g _ { i } ( \lambda ) & = f _ { i } ( \lambda ) + \hat { \mathfrak T } \theta _ { j } ( \mu ) w _ { j } ( \lambda ) g _ { j } ( \lambda ) & & ( \text{A.10} )$$

$$w _ { i } ( \lambda ) g _ { i } ( \lambda ) = w _ { i } ( \lambda ) b _ { i } ( \lambda ) + w _ { i } ( \lambda ) \hat { T } m _ { j } ( \mu ) g _ { j } ( \mu ).$$

Combining both equations we find:

$$w _ { i } ( \lambda ) = \frac { 1 } { b _ { i } ( \lambda ) } \left ( f _ { i } ( \lambda ) + \hat { T } \theta _ { j } ( \mu ) g _ { j } ( \lambda ) w _ { j } ( \lambda ) - w _ { i } ( \lambda ) \hat { T } m _ { j } ( \mu ) g _ { j } ( \mu ) \right ). \quad \text{(A.12)}$$

This is a new fixed point equation for wi ( λ ) . It is motivated by the equation for the effective velocity, which is obtained in the special case f i ( λ ) = E i ′ ( λ ) , b i ( λ ) = P ′ i ( λ ) and mi ( λ ) = ni ( λ ) . Intuitively speaking we need that ζ i ( λ ) decays fast for large i , but at the same time b i ( λ ) should decay slowly, such that (A.12) converges.

For the gapped phase we used g i ( λ ) = 1 2 / i and mi ( λ ) = 1, which happened to make (A.12) convergent with initial guess f i ( λ / ) b i ( λ ) .

Similarly, one can derive an equation for ˜ Ki ( t , x , λ ) = Ki ( t , x , λ / ) g i ( λ ) :

$$\tilde { K } _ { i } ( t, x, \lambda ) = \frac { 1 } { b _ { i } ( \lambda ) } \left ( P _ { i } ^ { \prime } ( \lambda ) x - E _ { i } ^ { \prime } ( \lambda ) t - \tilde { K } _ { i } ( t, x, \lambda ) \hat { T } m _ { j } ( \mu ) g _ { j } ( \mu ) + \hat { T } [ \hat { N } _ { 0 ; j } ( g _ { j } ( \mu ) \tilde { K } _ { j } ( t, x, \mu ), \mu ) ] \right ),$$

which for the partitioning protocol becomes, for g i ( λ &gt; ) 0:

$$\tilde { K } _ { i } ( t, x, \lambda ) & = \frac { 1 } { b _ { i } ( \lambda ) } \Big ( P _ { i } ^ { \prime } ( \lambda ) x - E _ { i } ^ { \prime } ( \lambda ) t - \tilde { K } _ { i } ( t, x, \lambda ) \hat { T } m _ { j } ( \mu ) g _ { j } ( \mu ) \\ & + \hat { T } [ \theta _ { \text{sgn} ( \tilde { K } _ { j } ( t, x, \mu ) ) ; i j } ( \mu ) g _ { j } ( \mu ) \tilde { K } _ { j } ( t, x, \mu ) ] \Big ).$$

This equation, with the above stated g i ( λ ) = 1 2 / i and mi ( λ ) = 1, is solved by iteration, starting from P ′ dr sgn ( x ) , i ( λ ) b i ( λ ) x , to obtain the solution to the GHD equation in the gapped case.

## B Derivation of the GGE rapidity distribution functions

In this section we describe the derivation of the GGE rapidity distribution functions for the quasiparticles and holes ρ n ( λ ) , ρ h n ( λ ) , which emerges after local relaxation from the initial states involved in our bipartition protocol.

The derivation follows from the boundary quantum transfer matrix (BQTM) approach explained in Ref. [ 86 ] for the case of a continuous Hamiltonian evolution, see also [ 87,88,114116 , and extended in ] [ 33 ] to the circuit dynamics.

## B.1 The boundary quantum transfer matrix construction

Our starting point is a continuous family of initial states

$$| \Psi _ { 0 } ( z ) \rangle = & \left ( \frac { z ^ { 1 / 2 } | 0 1 \rangle - z ^ { - 1 / 2 } | 1 0 \rangle } { \sqrt { | z | + | z | ^ { - 1 } } } \right ) ^ { \otimes L / 2 }$$

which interpolates between the various initial states considered in our partitioning protocol. More concretely, Néel, anti-Néel and dimer states are recovered for z = 0, z = ∞ and z = 1 respectively.

The BQTM approach allows to compute the GGE rapidity distributions from initial states of the form (B.1) by recasting the discrete time evolution in terms of auxiliary transfer matrices acting in the space direction, and where the initial state enters as a boundary condition.

We refer to the previous works [ 33,86-88,114 ] for a complete definition of the boundary QTM. Its boundaries are specified by a set of parameters, which in the present case are given by

$$\xi _ { \pm } = \arctan \left ( \frac { z + 1 } { z - 1 } \tan \left ( \frac { x } { 2 } \pm \frac { \gamma } { 2 } \right ) \right )$$

(as oppposed to the case of most generic two-site product states considered in [ 87 , the states ] (B.1) correspond to diagonal boundary conditions, so all the other parameters vanish).

In the case of Hamiltonian evolution, the BQTM was originally introduced to compute a quantity called the Loschmidt echo [ 87, 114 , resulting in a transfer matrix acting on ] N auxiliary spins (where N is a introduced as a Trotterization number for the time evolution). Physical quantities of interest were then related to the dominant eigenvalue of this auxiliary transfer matrix, expressed in terms of a set of auxiliary Bethe roots. In particular, the GGE rapidity distributions are obtained by taking a 'zero time limit' in the BQTM generating the Loschmidt echo. In this limit, the auxiliary Bethe roots collapse to some known analytical value and their contribution to the dominant eigenvalue cancels. As a result, the GGE rapidity distributions can then be simply recovered from a BQTM with N = 0 sites, namely a scalar. The same observations hold in the circuit case (where a BQTM with N = 2 was used in our previous work [ 33 ), and in this work we shall therefore work from the start with the ] N = 0 version of the BQTM, which reads

$$\mathcal { T } ( u ) = \omega _ { 1 } ( u ) + \omega _ { 2 } ( u ) \,,$$

where

$$\omega _ { 1 } ( u ) \ & = \ \frac { \sin ( 2 u + \gamma ) \sin ( u + \xi ^ { + } - \frac { \gamma } { 2 } ) \sin ( u + \xi ^ { - } - \frac { \gamma } { 2 } ) } { \sin ( 2 u ) } \,, \\ \omega _ { 2 } ( u ) \ & = \ \frac { \sin ( 2 u - \gamma ) \sin ( u - \xi ^ { + } + \frac { \gamma } { 2 } ) \sin ( u - \xi ^ { - } + \frac { \gamma } { 2 } ) } { \sin ( 2 u ) } \,.$$

## B.2 The T-system and Y-system

The (scalar) transfer matrices T ( u ) are part of a hierarchy of commuting transfer matrices T i ( u ) , satisfying the recursive relations [ 87 ]

$$\mathcal { T } _ { 0 } \ & = \ 1, \\ \mathcal { T } _ { 1 } \ & = \ \mathcal { T }, \\ \mathcal { T } _ { j } \ & = \ \mathcal { T } _ { j - 1 } ^ { - } \mathcal { T } _ { 1 } ^ { [ j - 1 ] } - f ^ { [ j - 3 ] } \mathcal { T } _ { j - 2 } ^ { [ - 2 ] }, \quad j \geq 2,$$

where for any function F u ( ) wehaveintroduced the short-hand notation F [ k ] ( u ) ≡ F u ( + ik γ/ 2 , ) and where f ( u -γ/ 2 ) = ω 1 ( u + γ/ 2 ) ω 2 ( u -γ/ 2 . In the more general case where the transfer ) matrix acts on a chain of N spins, or for more general (e.g. non-diagonal) boundary conditions, the definition of f is changed, but the general structure (B.5) remains the same [ 87 . ]

The transfer matrices T j satisfy a set of relations known as the T-system [ 87, 117, 118 , ] and which read

$$\mathcal { T } _ { j } ^ { [ m ] } \mathcal { T } _ { j } ^ { [ - m ] } = \mathcal { T } _ { j + m } \mathcal { T } _ { j - m } + \Phi _ { j - m + 1 } \mathcal { T } _ { m - 1 } ^ { [ j + 1 ] } \mathcal { T } _ { m - 1 } ^ { [ - j - 1 ] } \,,$$

where we have introduced the function Φ j = Q j f [ 2 k - -2 j ]

.

One then introduces the family of Y-functions Y 0 = 0, and l = 1

$$Y _ { j } ( \lambda ) = \frac { \mathcal { T } _ { j - 1 } ( u ) \mathcal { T } _ { j + 1 } ( u ) } { \Phi _ { j } ( u ) } \,, \quad j \geq 1 \,.$$

As a consequence of the T-system, the Y-functions obey the following set of relations, known as Y-system [ 117 ]

$$Y _ { j } \left ( u + \frac { \gamma } { 2 } \right ) Y _ { j } \left ( u - \frac { \gamma } { 2 } \right ) = \left [ 1 + Y _ { j + 1 } \left ( u \right ) \right ] \left [ 1 + Y _ { j - 1 } \left ( u \right ) \right ] \,,$$

and can therefore be obtained recursively from the knowledge of Y 1 ( λ ) , which is

$$1 + Y _ { 1 } ( u ) = \left ( 1 + \mathfrak { a } ( u - \gamma / 2 ) \right ) \left ( 1 + 1 / \mathfrak { a } ( u + \gamma / 2 ) \right ) \,,$$

where

$$\mathfrak { a } ( u ) = \frac { \omega _ { 1 } ( u ) } { \omega _ { 2 } ( u ) } = \frac { \sin ( 2 u + \gamma ) } { \sin ( 2 u - \gamma ) } \frac { \sin ( u + \xi _ { + } - \gamma / 2 ) \sin ( u + \xi _ { - } - \gamma / 2 ) } { \sin ( u - \xi _ { + } + \gamma / 2 ) \sin ( u - \xi _ { - } + \gamma / 2 ) } \,. \quad \quad ( B. 1 0 ) \\.. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \. \.$$

Note that for the dimer state, z = 1, the dependence in x vanishes (which can already be noted from the fact that ξ ± do not depend on x ).

As shown in previous works, the GGE functions η j ( λ ) are then obtained from the Yfunctions [ 87, 88 . ] For generic values of the anisotropy ∆ , in particular in the gapped phase | ∆ | &gt; 1 (that is, γ = i η η , ∈ Ρ ), the GGE is characterized by an infinite family { η j ( λ ) } j ≥ 1 , which are simply identified one by one with the Y functions, i.e. η j ( λ ) = Yj ( λ ) .

Things are however more subtle in the gapless phase | ∆ | &lt; 1, which is spanned by the socalled root-of-unity points γ/π ∈ Θ : for such values of γ the TBA truncates to a finite number of rapidity distribution functions { η j ( λ ) } j = 1,... Nb . Accordingly, the T- and Y- system can also be truncated to a finite number of equations, and we now turn to a description of the precise correspondence between the two, limiting ourselves, for simplicity, to values of the anisotropy of the form γ = π p + 1 , with p an integer ≥ 1.

## B.3 Truncation of the T- and Y- system at γ = π p + 1

At γ = π p + 1 there exists an additional linear relation between the transfer matrices. It generally takes the form [ 88 ]

$$\mathcal { T } _ { p + 1 } ( u ) = \alpha ( u ) \mathcal { T } _ { p - 1 } ( u ) + \beta ( u ) \,,$$

where, in the present case, the functions α ( u ) and β ( u ) are given by

$$\alpha ( u ) \ = \ \frac { \left ( \Omega ^ { [ - p - 2 ] } \Omega _ { 1 } ^ { [ - p - 2 ] } \Omega _ { 2 } ^ { [ - p - 2 ] } \right ) ^ { 2 } } { \Phi _ { p } } \,, \\ \quad \ \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Phi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi_ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ {p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p } \, \Psi _ { p }$$

$$( B. 1 2 )$$

$$\beta ( u ) \ = \ \Omega ^ { - [ p + 2 ] } \left ( ( \Omega _ { 1 } ^ { [ - p - 2 ] } ) ^ { 2 } + ( \Omega _ { 2 } ^ { [ - p - 2 ] } ) ^ { 2 } \right ) \,,$$

and where we have defined

$$\Omega ( u ) \ = \ \cot ( ( p + 1 ) u ) \cot ( ( p + 1 ) ( u + \pi / 2 ) ) \,, \quad.$$

cos

(

p

+

1

)(

u

±

GLYPH&lt;0&gt;

GLYPH&lt;1&gt;

GLYPH&lt;0&gt;

±

u

(

))

=

ξ

+ )

cos

(

p

+

1

)(

u

2 2

p

2

Ω

(

1,2

ξ

-

GLYPH&lt;1&gt;

)

.

(B.14)

As a result of (B.11), the Y-system can also be truncated: introducing the function

$$K ( u ) = \frac { \Omega ^ { [ - p - 2 ] } \Omega _ { 1 } ^ { [ - p - 2 ] } \Omega _ { 2 } ^ { [ - p - 2 ] } } { \Phi _ { p } } \mathcal { T } _ { p - 1 } ( u ) \,,$$

we have

$$Y _ { j } ^ { + } Y _ { j } ^ { - } & = ( 1 + Y _ { j + 1 } ) ( 1 + Y _ { j - 1 } ) \,, \quad j = 1, \dots p - 1 \,, \\ 1 + Y _ { p - 1 } & = K ^ { + } K ^ { - } \,, \quad \\ \quad \cdot \quad \cdot \quad \Omega _ { \cdot } ^ { [ - p - 2 ] }$$

$$1 + Y _ { p } = 1 + ( \mathfrak { b } + \mathfrak { b } ^ { - 1 } ) K + K ^ { 2 } \,, \quad \mathfrak { b } = \frac { \Omega _ { 1 } ^ { [ - p - 2 ] } } { \Omega _ { 2 } ^ { [ - p - 2 ] } } \,.$$

## B.4 The GGE densities at γ = π/ ( p + 1 )

Using classic techniques, the truncated Y-system can be turned into a set of non-linear integral equations (NLIE) [ 88 . Those can in turn be compared with the generic structure of the NLIE ] obeyed by the functions η j involved in the TBA solution, which can be found in Section 9.2 of [ 66 ] for the gapless XXZ Hamiltonian at finite temperature, and retain the same general structure for Generalized Gibbs Ensembles or for the discrete dynamics (only the explicit form of source terms gets modified). In order to match the form of the NLIE inherited from the Y-system with those of the TBA, we are led to the following identification

$$\eta _ { j } ( \lambda ) & = Y _ { j } ( i \lambda ) \,, \quad j = 1, \dots p - 1 \\ \eta _ { p } ( \lambda ) & = \mathfrak { b } ( i \lambda ) K ( i \lambda ) \\ \eta _ { p + 1 } ( \lambda ) & = \mathfrak { b } ( i \lambda ) / K ( i \lambda ) \\ \intertext { m a t i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i n d i } \,,$$

For concreteness, we give in the main text the explicit expressions of the functions η j ( λ ) for p = 2 and p = 3, see Eqs. (52) and (51).

## C Details on the free-fermion points of the model

In this Appendix we provide details on the techniques we have used to study the non-interacting points of the model, namely those in which the circuit can be mapped onto the Floquet dynamics of free fermions.

From Eqs. (75) and (76), we see that the unitary operators Ue ( τ ) , Uo ( τ ) are written as the exponential of quadratic Hamiltonians. This makes it possible to diagonalise the Floquet operator U ( τ ) by means of an appropriate Bogoliubov rotation [ 33 . In addition, the mapping ] to free fermions makes it possible to simulate efficiently the dynamics of the system, using the formalism of Gaussian states [ 97 . ] We used this method to produce the numerical data presented in Sec. 4.1, as we explain in the following.

We recall that Gaussian states are defined by the fact that they satisfy Wick's theorem [ 97 . ] Accordingly, given a pure Gaussian state | ψ 〉 with fixed particle number, it is fully determined by its covariance matrix

$$\Gamma _ { i j } = \langle \psi | c _ { i } ^ { \dagger } c _ { j } | \psi \rangle \,.$$

This observation drastically simplifies the problem, as Γ is a matrix whose linear dimensions grow linearly, not exponentially, with the system size.

The unitary updates of the system state | ψ 2 t + 1 〉 = Uo ( τ ) | ψ 2 t 〉 and | ψ 2 t 〉 = Ue ( τ ) | ψ 2 t -1 〉 correspond to a change in the covariance matrix which can be written as [ 97 ]

$$\Gamma ( 2 t + 1 ) = V _ { o } ^ { * } \Gamma ( 2 t ) V _ { o } ^ { T } \,,$$

$$\Gamma ( 2 t ) = V _ { e } ^ { * } \Gamma ( 2 t - 1 ) V _ { e } ^ { T } \,,$$

where ( ) · ∗ and ( ) · T denote complex conjugation and transposition, respectively. Here,

$$V _ { e } = e ^ { A _ { e } } \,, \quad V _ { o } = e ^ { A _ { o } } \,,$$

where Ae and Ao are L × L matrices whose non-zero elements are defined by

$$($$

$$( A _ { o } ) _ { 2 n, 2 n + 1 } = ( A _ { o } ) _ { 2 n + 1, 2 n } = - i \frac { \tau } { 2 } \,, & & ( \text{C.5} ) \\ ( A _ { e } ) _ { 2 n - 1, 2 n } = ( A _ { e } ) _ { 2 n, 2 n - 1 } = - i \frac { \tau } { 2 } \,. & & ( \text{C.6} ) \\ \cdot \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \cdot \, \cdot \, \cdot \, \cdot \, - \, \cdot \, \dots \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot$$

Wecan also map the initial states (42)-(44) into fermionic Fock states, which are Gaussian. This is true also for bipartite initial states, where the two halves of the system are prepared in different product states chosen among (42)-(44). In fact, it is straightforward to compute their covariance matrix. Given the initial matrix Γ , the evolution is computed numerically via iteration of Eqs. (C.2).

## References

- [ 1 ] H. F . Trotter, On the product of semi-groups of operators , Proc. American Math. Soc. 10 (4), 545 (1959).
- [ 2 ] M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics , J. Math. Phys. 32 (2), 400 (1991), doi:10.1063 / 1.529425.
- [ 3 ] I. M. Georgescu, S. Ashhab and F . Nori, Quantum simulation , Rev. Mod. Phys. 86 , 153 (2014), doi:10.1103 / RevModPhys.86.153.
- [ 4 ] S. Lloyd, Universal quantum simulators , Science 273 (5278), 1073 (1996), doi:10.1126 / science.273.5278.1073.
- [ 5 ] J. Preskill, Quantum computing in the nisq era and beyond , Quantum 2 , 79 (2018), doi:10.22331 / q-2018-08-06-79.
- [ 6 ] J. I. Cirac and P . Zoller, Goals and opportunities in quantum simulation , Nature Phys. 8 (4), 264 (2012), doi:10.1038 / nphys2275.
- [ 7 ] P . Arrighi, An overview of quantum cellular automata , Nat. Comp. 18 (4), 885 (2019), doi:10.1007 / s11047-019-09762-6.
- [ 8 ] T. Farrelly, A review of quantum cellular automata , Quantum 4 , 368 (2020), doi:10.22331 / q-2020-11-30-368.
- [ 9 ] I. Bloch, J. Dalibard and W . Zwerger, Many-body physics with ultracold gases , Rev. Mod. Phys. 80 , 885 (2008), doi:10.1103 / RevModPhys.80.885.
- [ 10 ] R. Blatt and C. F . Roos, Quantum simulations with trapped ions , Nature Physics 8 (4), 277 (2012), doi:10.1038 / nphys2252.
- [ 11 ] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms , Nature Phys. 16 (2), 132 (2020), doi:10.1038 / s41567-019-0733-z.
- [ 12 ] C. Monroe, W . C. Campbell, L.-M. Duan, Z.-X. Gong, A. V . Gorshkov, P . W . Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P . Richerme, C. Senko et al. , Programmable quantum simulations of spin systems with trapped ions , Rev. Mod. Phys. 93 , 025001 (2021), doi:10.1103 / RevModPhys.93.025001.

- [ 13 ] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer and P . Zoller, Practical quantum advantage in quantum simulation , Nature 607 (7920), 667 (2022), doi:10.1038 / s41586-022-04940-6.
- [ 14 ] M. Heyl, P . Hauke and P . Zoller, Quantum localization bounds trotter errors in digital quantum simulation , Science Adv. 5 (4), eaau8342 (2019), doi:10.1126 / sciadv.aau8342.
- [ 15 ] D. A. Abanin, W. De Roeck and F . Huveneers, Exponentially slow heating in periodically driven many-body systems , Phys. Rev. Lett. 115 , 256803 (2015), doi:10.1103 / PhysRevLett.115.256803.
- [ 16 ] D. A. Abanin, W . De Roeck, W . W . Ho and F . Huveneers, Effective hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems , Phys. Rev. B 95 , 014112 (2017), doi:10.1103 / PhysRevB.95.014112.
- [ 17 ] D. Abanin, W . De Roeck, W . W . Ho and F . Huveneers, A rigorous theory of many-body prethermalization for periodically driven and closed quantum systems , Comm. Math. Phys. 354 (3), 809 (2017), doi:10.1007 / s00220-017-2930-x.
- [ 18 ] T. Mori, T. Kuwahara and K. Saito, Rigorous bound on energy absorption and generic relaxation in periodically driven quantum systems , Phys. Rev. Lett. 116 , 120401 (2016), doi:10.1103 / PhysRevLett.116.120401.
- [ 19 ] T. Kuwahara, T. Mori and K. Saito, Floquet-magnus theory and generic transient dynamics in periodically driven many-body quantum systems , Ann. Phys. 367 , 96 (2016), doi:https: // doi.org / 10.1016 j.aop.2016.01.012. /
- [ 20 ] T. Ishii, T . Kuwahara, T . Mori and N. Hatano, Heating in integrable time-periodic systems , Phys. Rev. Lett. 120 , 220602 (2018), doi:10.1103 / PhysRevLett.120.220602.
- [ 21 ] L. M. Sieberer, T . Olsacher, A. Elben, M. Heyl, P . Hauke, F . Haake and P . Zoller, Digital quantum simulation, trotter errors, and quantum chaos of the kicked top , npj Quantum Inf. 5 (1), 1 (2019), doi:10.1038 / s41534-019-0192-5.
- [ 22 ] C. Kargi, J. P . Dehollain, F . Henriques, L. M. Sieberer, T . Olsacher, P . Hauke, M. Heyl, P . Zoller and N. K. Langford, Quantum chaos and universal trotterisation behaviours in digital quantum simulations , arXiv:2110.11113 (2021).
- [ 23 ] K. Chinni, M. H. Muñoz Arias, I. H. Deutsch and P . M. Poggi, Trotter errors from dynamical structural instabilities of floquet maps in quantum simulation , PRX Quantum 3 , 010351 (2022), doi:10.1103 / PRXQuantum.3.010351.
- [ 24 ] M. Rigol, V . Dunjko and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems , Nature 452 (7189), 854 (2008), doi:10.1038 / nature06838.
- [ 25 ] L. D'Alessio, Y. Kafri, A. Polkovnikov and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics , Adv. Phys. 65 (3), 239 (2016), doi:10.1080 / 00018732.2016.1198134.
- [ 26 ] M. Rigol, V . Dunjko, V . Yurovsky and M. Olshanii, Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons , Phys. Rev. Lett. 98 , 050405 (2007), doi:10.1103 / PhysRevLett.98.050405.

- [ 27 ] L. Vidmar and M. Rigol, Generalized gibbs ensemble in integrable lattice models , J. Stat. Mech. 2016 (6), 064007 (2016).
- [ 28 ] F . H. Essler and M. Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains , J. Stat. Mech. 2016 (6), 064002 (2016).
- [ 29 ] M. Vanicat, L. Zadnik and T. Prosen, Integrable trotterization: Local conservation laws and boundary driving , Phys. Rev. Lett. 121 , 030606 (2018), doi:10.1103 / PhysRevLett.121.030606.
- [ 30 ] M. Ljubotina, L. Zadnik and T. Prosen, Ballistic spin transport in a periodically driven integrable quantum system , Phys. Rev. Lett. 122 , 150605 (2019), doi:10.1103 / PhysRevLett.122.150605.
- [ 31 ] M. Medenjak, T. Prosen and L. Zadnik, Rigorous bounds on dynamical response functions and time-translation symmetry breaking , SciPost Phys. 9 , 003 (2020), doi:10.21468 / SciPostPhys.9.1.003.
- [ 32 ] G. Giudice, G. Giudici, M. Sonner, J. Thoenniss, A. Lerose, D. A. Abanin and L. Piroli, Temporal entanglement, quasiparticles, and the role of interactions , Phys. Rev. Lett. 128 , 220401 (2022), doi:10.1103 / PhysRevLett.128.220401.
- [ 33 ] E. Vernier, B. Bertini, G. Giudici and L. Piroli, Integrable digital quantum simulation: Generalized gibbs ensembles and trotter transitions , Phys. Rev. Lett. 130 , 260401 (2023), doi:10.1103 / PhysRevLett.130.260401.
- [ 34 ] E. Vernier, H.-C. Yeh, L. Piroli and A. Mitra, Strong zero modes in integrable quantum circuits , arXiv:2401.12305 (2024).
- [ 35 ] L. Zadnik, M. Ljubotina, Ž. Krajnik, E. Ilievski and T . Prosen, Quantum many-body spin ratchets , arXiv:2406.01571 (2024).
- [ 36 ] I. L. Aleiner, Bethe ansatz solutions for certain periodic quantum circuits , Ann. Phys. 433 , 168593 (2021), doi:10.1016 / j.aop.2021.168593.
- [ 37 ] Y. Miao, V . Gritsev and D. V . Kurlov, The floquet baxterisation , arXiv:2206.15142 (2022).
- [ 38 ] P . W. Claeys, J. Herzog-Arbeitman and A. Lamacraft, Correlations and commuting transfer matrices in integrable unitary circuits , SciPost Phys. 12 , 007 (2022), doi:10.21468 / SciPostPhys.12.1.007.
- [ 39 ] C. Paletta and T. Prosen, Integrability of open boundary driven quantum circuits , arXiv:2406.12695 (2024).
- [ 40 ] A. Morvan, T. I. Andersen, X. Mi, C. Neill, A. Petukhov, K. Kechedzhi, D. Abanin, R. Acharya, F . Arute, K. Arya et al. , Formation of robust bound states of interacting photons , arXiv:2206.05254 (2022).
- [ 41 ] 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: Math. Theor. 56 (16), 165301 (2023), doi:10.1088 / 1751-8121 acc369. /
- [ 42 ] N. Keenan, N. Robertson, T. Murphy, S. Zhuk and J. Goold, Evidence of kardar-parisizhang scaling on a digital quantum simulator , arXiv:2208.12243 (2022).
- [ 43 ] V . E. Korepin, N. M. Bogoliubov and A. G. Izergin, Quantum inverse scattering method and correlation functions , vol. 3, Cambridge university press (1997).

- [ 44 ] T. Antal, Z. Rácz, A. Rákos and G. M. Schütz, Transport in the XX chain at zero temperature: Emergence of flat magnetization profiles , Phys. Rev. E 59 , 4912 (1999), doi:10.1103 / PhysRevE.59.4912.
- [ 45 ] W. H. Aschbacher and C.-A. Pillet, Non-equilibrium steady states of the xy chain , J. Stat. Phys 112 , 1153-1175 (2003), doi:10.1023 / A:1024619726273.
- [ 46 ] W. H. Aschbacher and J.-M. Barbaroux, Out of equilibrium correlations in the xy chain , Lett. Math. Phys. 77 , 11 (2006), doi:10.1007 / s11005-006-0049-7.
- [ 47 ] V . Alba, B. Bertini, M. Fagotti, L. Piroli and P . Ruggiero, Generalized-hydrodynamic approach to inhomogeneous quenches: correlations, entanglement and quantum effects , J. Stat. Mech. 2021 (11), 114004 (2021), doi:10.1088 / 1742-5468 ac257d. /
- [ 48 ] P . Calabrese and J. Cardy, Time dependence of correlation functions following a quantum quench , Phys. Rev. Lett. 96 , 136801 (2006), doi:10.1103 / PhysRevLett.96.136801.
- [ 49 ] B. Bertini, M. Collura, J. De Nardis and M. Fagotti, Transport in out-of-equilibrium x xz chains: Exact profiles of charges and currents , Phys. Rev. Lett. 117 , 207201 (2016), doi:10.1103 / PhysRevLett.117.207201.
- [ 50 ] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium , Phys. Rev. X 6 , 041065 (2016), doi:10.1103 / PhysRevX.6.041065.
- [ 51 ] B. Doyon and H. Spohn, Dynamics of hard rods with initial domain wall state , J. Stat. Mech. 2017 (7), 073210 (2017), doi:10.1088 / 1742-5468 aa7abf. /
- [ 52 ] A. De Luca, M. Collura and J. De Nardis, Nonequilibrium spin transport in integrable spin chains: Persistent currents and emergence of magnetic domains , Phys. Rev. B 96 , 020403 (2017), doi:10.1103 / PhysRevB.96.020403.
- [ 53 ] L. Piroli, J. De Nardis, M. Collura, B. Bertini and M. Fagotti, Transport in out-ofequilibrium xxz chains: Nonballistic behavior and correlation functions , Phys. Rev. B 96 , 115124 (2017), doi:10.1103 / PhysRevB.96.115124.
- [ 54 ] M. Mestyán, B. Bertini, L. Piroli and P . Calabrese, Spin-charge separation effects in the low-temperature transport of one-dimensional fermi gases , Phys. Rev. B 99 , 014305 (2019), doi:10.1103 / PhysRevB.99.014305.
- [ 55 ] V . Alba, B. Bertini and M. Fagotti, Entanglement evolution and generalised hydrodynamics: interacting integrable systems , SciPost Phys. 7 , 5 (2019), doi:10.21468 / SciPostPhys.7.1.005.
- [ 56 ] B. Bertini, L. Piroli and P . Calabrese, Universal broadening of the light cone in low-temperature transport , Phys. Rev. Lett. 120 , 176801 (2018), doi:10.1103 / PhysRevLett.120.176801.
- [ 57 ] B. Bertini and L. Piroli, Low-temperature transport in out-of-equilibrium XXZ chains , J. Stat. Mech. 2018 , 033104 (2018), doi:10.1088 / 1742-5468 aab04b. /
- [ 58 ] A. Bastianello, B. Doyon, G. Watts and T. Yoshimura, Generalized hydrodynamics of classical integrable field theory: the sinh-gordon model , SciPost Phys. 4 , 45 (2018), doi:10.21468 / SciPostPhys.4.6.045.

- [ 59 ] B. Doyon, Generalized hydrodynamics of the classical toda system , J. Math. Phys. 60 , 073302 (2019).
- [ 60 ] V . B. Bulchandani, X. Cao and J. E. Moore, Kinetic theory of quantum and classical toda lattices , J. Phys. A: Math. Theor. 52 , 33LT01 (2019), doi:10.1088 / 1751-8121 ab2cf0. /
- [ 61 ] Y. Nozawa and H. Tsunetsugu, Generalized hydrodynamic approach to charge and energy currents in the one-dimensional hubbard model , Phys. Rev. B 101 , 035121 (2020), doi:10.1103 / PhysRevB.101.035121.
- [ 62 ] Y. Nozawa and H. Tsunetsugu, Generalized hydrodynamics study of the one-dimensional hubbard model: Stationary clogging and proportionality of spin, charge, and energy currents , Phys. Rev. B 103 , 035130 (2021), doi:10.1103 / PhysRevB.103.035130.
- [ 63 ] M. Mestyàn and V . Alba, Molecular dynamics simulation of entanglement spreading in generalized hydrodynamics , SciPost Phys. 8 , 55 (2020), doi:10.21468 / SciPostPhys.8.4.055.
- [ 64 ] S. Wang, X. Yin, Y.-Y. Chen, Y. Zhang and X.-W . Guan, Emergent ballistic transport of bose-fermi mixtures in one dimension , J. Phys. A: Math. Theor. 53 (46), 464002 (2020), doi:10.1088 / 1751-8121 abc128. /
- [ 65 ] 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.
- [ 66 ] M. Takahashi, Thermodynamics of one-dimensional solvable models , Cambridge University Press (2005).
- [ 67 ] M. Mestyán and B. Pozsgay, Short distance correlators in the xxz spin chain for arbitrary string distributions , J. Stat. Mech. 2014 (9), P09020 (2014), doi:10.1088 / 17425468 2014 09 P09020. / / /
- [ 68 ] L. Piroli, E. Vernier and P . Calabrese, Exact steady states for quantum quenches in integrable heisenberg spin chains , Phys. Rev. B 94 , 054313 (2016), doi:10.1103 / PhysRevB.94.054313.
- [ 69 ] L. Piroli, E. Vernier, P . Calabrese and M. Rigol, Correlations and diagonal entropy after quantum quenches in xxz chains , Phys. Rev. B 95 , 054308 (2017), doi:10.1103 / PhysRevB.95.054308.
- [ 70 ] A. Bastianello, L. Piroli and P . Calabrese, Exact local correlations and full counting statistics for arbitrary states of the one-dimensional interacting bose gas , Phys. Rev. Lett. 120 , 190601 (2018), doi:10.1103 / PhysRevLett.120.190601.
- [ 71 ] B. Pozsgay, Algebraic construction of current operators in integrable spin chains , Phys. Rev. Lett. 125 , 070602 (2020), doi:10.1103 / PhysRevLett.125.070602.
- [ 72 ] M. Borsi, B. Pozsgay and L. Pristyák, Current operators in integrable models: a review , J. Stat. Mech. 2021 (9), 094001 (2021), doi:10.1088 / 1742-5468 ac0f6b. /
- [ 73 ] L. Pristyák and B. Pozsgay, Current mean values in the XYZ model , SciPost Phys. 14 , 158 (2023), doi:10.21468 / SciPostPhys.14.6.158.
- [ 74 ] T. Platini and D. Karevski, Relaxation in the xx quantum chain , J. Phys. A: Math. Theor. 40 (8), 1711 (2007), doi:10.1088 / 1751-8113 40 8 002. / / /

- [ 75 ] J. Lancaster and A. Mitra, Quantum quenches in an x xz spin chain from a spatially inhomogeneous initial state , Phys. Rev. E 81 , 061134 (2010), doi:10.1103 / PhysRevE.81.061134.
- [ 76 ] V . Eisler and Z. Rácz, Full counting statistics in a propagating quantum front and random matrix spectra , Phys. Rev. Lett. 110 , 060602 (2013), doi:10.1103 / PhysRevLett.110.060602.
- [ 77 ] A. De Luca, J. Viti, D. Bernard and B. Doyon, Nonequilibrium thermal transport in the quantum ising chain , Phys. Rev. B 88 , 134301 (2013), doi:10.1103 / PhysRevB.88.134301.
- [ 78 ] B. Doyon, A. Lucas, K. Schalm and M. Bhaseen, Non-equilibrium steady states in the klein-gordon theory , J. Phys. A: Math. Theor. 48 (9), 095002 (2015), doi:10.1088 / 1751-8113 48 9 095002. / / /
- [ 79 ] N. Allegra, J. Dubail, J.-M. Stéphan and J. Viti, Inhomogeneous field theory inside the arctic circle , J. Stat. Mech. 2016 (5), 053108 (2016), doi:10.1088 / 17425468 2016 05 053108. / / /
- [ 80 ] D. Bernard and B. Doyon, Non-equilibrium steady states in conformal field theory , In Ann. Henri Poincaré , vol. 16, pp. 113-161. Springer, doi:10.1007 / s00023-014-0314-8 (2015).
- [ 81 ] D. Bernard and B. Doyon, A hydrodynamic approach to non-equilibrium conformal field theories , J. Stat. Mech. 2016 (3), 033104 (2016), doi:10.1088 / 17425468 2016 03 033104. / / /
- [ 82 ] J. Dubail, J.-M. Stéphan, J. Viti and P . Calabrese, Conformal field theory for inhomogeneous one-dimensional quantum systems: the example of non-interacting fermi gases , SciPost Phys. 2 , 002 (2017), doi:10.21468 / SciPostPhys.2.1.002.
- [ 83 ] G. Perfetto and A. Gambassi, Ballistic front dynamics after joining two semi-infinite quantum ising chains , Phys. Rev. E 96 , 012138 (2017), doi:10.1103 / PhysRevE.96.012138.
- [ 84 ] J. Dubail, J.-M. Stéphan and P . Calabrese, Emergence of curved light-cones in a class of inhomogeneous luttinger liquids , SciPost Phys. 3 , 019 (2017), doi:10.21468 / SciPostPhys.3.3.019.
- [ 85 ] L. Vidmar, D. Iyer and M. Rigol, Emergent eigenstate solution to quantum dynamics far from equilibrium , Phys. Rev. X 7 , 021012 (2017), doi:10.1103 / PhysRevX.7.021012.
- [ 86 ] L. Piroli, B. Pozsgay and E. Vernier, What is an integrable quench? , Nucl. Phys. B 925 , 362 (2017), doi:10.1016 / j.nuclphysb.2017.10.012.
- [ 87 ] L. Piroli, B. Pozsgay and E. Vernier, From the quantum transfer matrix to the quench action: the loschmidt echo in xxz heisenberg spin chains , J. Stat. Mech. 2017 (2), 023106 (2017), doi:10.1088 / 1742-5468 aa5d1e. /
- [ 88 ] L. Piroli, B. Pozsgay and E. Vernier, Non-analytic behavior of the loschmidt echo in xxz spin chains: Exact results , Nucl. Phys. B 933 , 454 (2018), doi:10.1016 / j.nuclphysb.2018.06.015.
- [ 89 ] B. Wouters, J. De Nardis, M. Brockmann, D. Fioretto, M. Rigol and J.-S. Caux, Quenching the anisotropic heisenberg chain: Exact solution and generalized gibbs ensemble predictions , Phys. Rev. Lett. 113 , 117202 (2014), doi:10.1103 / PhysRevLett.113.117202.

- [ 90 ] M. Brockmann, B. Wouters, D. Fioretto, J. De Nardis, R. Vlijm and J.-S. Caux, Quench action approach for releasing the néel state into the spin-1 / 2 xxz chain , J. Stat. Mech. 2014 (12), P12009 (2014), doi:10.1088 / 1742-5468 2014 12 P12009. / / /
- [ 91 ] B. Pozsgay, M. Mestyán, M. A. Werner, M. Kormos, G. Zaránd and G. Takács, Correlations after quantum quenches in the x xz spin chain: Failure of the generalized gibbs ensemble , Phys. Rev. Lett. 113 , 117203 (2014), doi:10.1103 / PhysRevLett.113.117203.
- [ 92 ] M. Mestyán, B. Pozsgay, G. Takács and M. Werner, Quenching the xxz spin chain: quench action approach versus generalized gibbs ensemble , J. Stat. Mech. 2015 (4), P04001 (2015), doi:10.1088 / 1742-5468 2014 12 P12009. / / /
- [ 93 ] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F . H. L. Essler and T . Prosen, Complete generalized gibbs ensembles in an interacting theory , Phys. Rev. Lett. 115 , 157201 (2015), doi:10.1103 / PhysRevLett.115.157201.
- [ 94 ] J.-S. Caux, The quench action , J. Stat. Mech. 2016 (6), 064006 (2016).
- [ 95 ] F . S. Møller and J. Schmiedmayer, Introducing ifluid: a numerical framework for solving hydrodynamical equations in integrable models , SciPost Phys. 8 , 041 (2020), doi:10.21468 / SciPostPhys.8.3.041.
- [ 96 ] F . Hübner and B. Doyon, A new quadrature for the generalized hydrodynamics equation and absence of shocks in the lieb-liniger model , arXiv:2406.18322 (2024).
- [ 97 ] S. Bravyi, Lagrangian representation for fermionic linear optics , arXiv quant-ph / 0404180 (2004).
- [ 98 ] M. Fagotti, M. Collura, F . H. L. Essler and P . Calabrese, Relaxation after quantum quenches in the spin1 2 heisenberg xxz chain , Phys. Rev. B 89 , 125101 (2014), doi:10.1103 / PhysRevB.89.125101.
- [ 99 ] B. Bertini and M. Fagotti, Determination of the nonequilibrium steady state emerging from a defect , Phys. Rev. Lett. 117 , 130402 (2016), doi:10.1103 / PhysRevLett.117.130402.
- [ 100 ] M. Fagotti, Global quenches after localized perturbations , Phys. Rev. Lett. 128 , 110602 (2022), doi:10.1103 / PhysRevLett.128.110602.
- [ 101 ] M. Fagotti, V . Mari´ c and L. Zadnik, Nonequilibrium symmetry-protected topological order: Emergence of semilocal gibbs ensembles , Phys. Rev. B 109 , 115117 (2024), doi:10.1103 / PhysRevB.109.115117.
- [ 102 ] V . Alba and P . Calabrese, Entanglement and thermodynamics after a quantum quench in integrable systems , PNAS 114 (30), 7947 (2017), doi:10.1073 / pnas.1703516114.
- [ 103 ] K. Temme, S. Bravyi and J. M. Gambetta, Error mitigation for short-depth quantum circuits , Phys. Rev. Lett. 119 , 180509 (2017), doi:10.1103 / PhysRevLett.119.180509.
- [ 104 ] S. McArdle, X. Yuan and S. Benjamin, Error-mitigated digital quantum simulation , Phys. Rev. Lett. 122 , 180501 (2019), doi:10.1103 / PhysRevLett.122.180501.
- [ 105 ] B. Koczor, Exponential error suppression for near-term quantum devices , Phys. Rev. X 11 , 031057 (2021), doi:10.1103 / PhysRevX.11.031057.
- [ 106 ] P . Czarnik, A. Arrasmith, P . J. Coles and L. Cincio, Error mitigation with clifford quantumcircuit data , Quantum 5 , 592 (2021), doi:10.22331 / q-2021-11-26-592.

- [ 107 ] A. Strikis, D. Qin, Y. Chen, S. C. Benjamin and Y. Li, Learning-based quantum error mitigation , PRX Quantum 2 , 040330 (2021), doi:10.1103 / PRXQuantum.2.040330.
- [ 108 ] C. Piveteau, D. Sutter and S. Woerner, Quasiprobability decompositions with reduced sampling overhead , npj Quantum Information 8 (1), 12 (2022), doi:10.1038 / s41534022-00517-3.
- [ 109 ] Y. Guo and S. Yang, Quantum error mitigation via matrix product operators , PRX Quantum 3 , 040313 (2022), doi:10.1103 / PRXQuantum.3.040313.
- [ 110 ] E. Van Den Berg, Z. K. Minev, A. Kandala and K. Temme, Probabilistic error cancellation with sparse pauli-lindblad models on noisy quantum processors , Nature Phys. 19 (8), 1116 (2023), doi:10.1038 / s41567-023-02042-2.
- [ 111 ] S. Filippov, M. Leahy, M. A. Rossi and G. García-Pérez, Scalable tensor-network error mitigation for near-term quantum computing , arXiv:2307.11740 (2023).
- [ 112 ] S. Mangini, M. Cattaneo, D. Cavalcanti, S. Filippov, M. A. Rossi and G. García-Pérez, Tensor network noise characterization for near-term quantum computers , arXiv:2402.08556 (2024).
- [ 113 ] King's College London. (2022). King's Computational Research, Engineering and Technology Environment (CREATE). Retrieved June 30, 2023, from https: // doi.org / 10.18742 rnvf-m076 / .
- [ 114 ] B. Pozsgay, The dynamical free energy and the loschmidt echo for a class of quantum quenches in the heisenberg spin chain , J. Stat. Mech. 2013 (10), P10028 (2013), doi:10.1088 / 1742-5468 2013 10 P10028. / / /
- [ 115 ] L. Piroli, E. Vernier, P . Calabrese and B. Pozsgay , Integrable quenches in nested spin chains i: the exact steady states , J. Stat. Mech. 2019 (6), 063103 (2019), doi:10.1088 / 17425468 ab1c51. /
- [ 116 ] L. Piroli, E. Vernier, P . Calabrese and B. Pozsgay, Integrable quenches in nested spin chains ii: fusion of boundary transfer matrices , J. Stat. Mech. 2019 (6), 063104 (2019), doi:10.1088 / 1742-5468 ab1c52. /
- [ 117 ] A. Kuniba, T. Nakanishi and J. Suzuki, T-systems and y-systems in integrable systems , J. Phys. A: Math. Theor. 44 (10), 103001 (2011), doi:10.1088 / 1742-5468 aa5d1e. /
- [ 118 ] Y.-k. Zhou, Row transfer matrix functional relations for baxter's eight-vertex and six-vertex models with open boundaries via more general reflection matrices , Nucl. Phys. B 458 (3), 504 (1996), doi:10.1016 / 0550-3213(95)00553-6.
Friedrich Hübner Eric Vernier Lorenzo Piroli
Generalized hydrodynamics of integrable quantum circuits
Quantum circuits make it possible to simulate the continuous-time dynamics of
a many-body Hamiltonian by implementing discrete Trotter steps of duration
$\tau$. However, when $\tau$ is sufficiently large, the discrete dynamics
exhibit qualitative differences compared to the original evolution, potentially
displaying novel features and many-body effects. We study an interesting
example of this phenomenon, by considering the integrable Trotterization of a
prototypical integrable model, the XXZ Heisenberg spin chain. We focus on the
well-known bipartition protocol, where two halves of a large system are
prepared in different macrostates and suddenly joined together, yielding
non-trivial nonequilibrium dynamics. Building upon recent results and adapting
the generalized hydrodynamics (GHD) of integrable models, we develop an exact
large-scale description of an explicit one-dimensional quantum-circuit setting,
where the input left and right qubits are initialized in two distinct product
states. We explore the phenomenology predicted by the GHD equations, which
depend on the Trotter step and the gate parameters. In some phases of the
parameter space, we show that the quantum-circuit large-scale dynamics is
qualitatively different compared to the continuous-time evolution. In
particular, we find that a single microscopic defect at the junction, such as
the addition of a single qubit, may change the nonequilibrium macrostate
appearing at late time.