The Quench Action

Upon quantization, the eigenstates of a Bethe Ansatz-solvable model [20–23] are uniquely labeled by a set of quantum numbers {I}. For simplicity, we consider a model in which only one type of particle is present (the generalization to cases where multiple types of string states are possible, e.g. for spin chains (see discussion around equation (115)), is completely straightforward, in which case the quantum numbers carry a 'particle type' index). If the model is in the continuum, the possible quantum numbers are unbounded; on the lattice, they obtain type-specific limits; we do not specify these model-dependent details at this stage. A resolution of the identity operator in Fock space can then be formally written as¹

$$\begin{eqnarray}&&\mathbf{1}=\underset{\left\{I\right\}}{\sum}|\left\{I\right\}\rangle \langle \left\{I\right\}|=\underset{N}{\sum}\underset{{{\left\{I\right\}}_{N}}}{\sum}|{{\left\{I\right\}}_{N}}\rangle \langle {{\left\{I\right\}}_{N}}|\end{eqnarray} \tag{ 3 }$$

in which the summation over N (the number of particles) runs over all allowed values. We use the bra-ket notation for eigenstates, namely $|{{\left\{I\right\}}_{N}}\rangle$ represents the normalized N-particle Bethe Ansatz state

$$\begin{eqnarray}&&|{{\left\{I\right\}}_{N}}\rangle \equiv \mathcal{N}\left({{\left\{I\right\}}_{N}}\right)\underset{j=1}{\overset{N}{\prod}}B\left({{\lambda}_{j}}\left({{\left\{I\right\}}_{N}}\right)\right)|0\rangle \end{eqnarray} \tag{ 4 }$$

obtained by acting on the reference state $|0\rangle$ with raising operators $B\left(\lambda \right)$ (see section 2.3) evaluated at the solutions of the Bethe equations in terms of rapidities ${{\lambda}_{j}}$ corresponding to quantum numbers {I}_N,

$$\begin{eqnarray}&&{{\theta}_{\text{kin}}}\left({{\lambda}_{j}}\right)+\frac{1}{L}\underset{l=1}{\overset{N}{\sum}}{{\theta}_{\text{scat}}}\left({{\lambda}_{j}}-{{\lambda}_{l}}\right)=\frac{2\pi}{L}{{I}_{j}},\quad j=1,...,N,\end{eqnarray} \tag{ 5 }$$

where L is the system size in dimensionless units and ${{\theta}_{\text{kin}}}$ and ${{\theta}_{\text{scat}}}$ are model-dependent kinetic and scattering kernels. In (4), $\mathcal{N}\left({{\left\{I\right\}}_{N}}\right)$ is the state normalization factor obtained from the appropriate Gaudin determinant. For definiteness, in (4), we have explicitly written each rapidity as being a function of the whole set of quantum numbers. This fact, which the reader should not forget, will be kept implicit in the following.

The general idea behind the thermodynamic limit is (as far as individual states are concerned) to replace sets of rapidities by distribution functions, and (and as far as state ensembles are concerned) to express the summations over eigenstates by a form of functional integration, with an eye towards the obtention of analytical results. The first step in this direction is to think of the quantum numbers scaled by the system size I_j/L in terms of a real variable x taking values on the real line. We associate to each eigenstate a root distribution ${{\rho}_{\text{dist}}}$ with support on the real axis

$$\begin{eqnarray}&&|\left\{I\right\}\rangle \Leftrightarrow {{\rho}_{\text{dist}}}\left(x;\left\{I\right\}\right),\end{eqnarray} \tag{ 6 }$$

this distribution in x space being defined as

$$\begin{eqnarray}&&{{\rho}_{\text{dist}}}\left(x;\left\{I\right\}\right)=\frac{1}{L}\underset{j}{\sum}\delta \left(x-{{x}_{j}}\right),\quad {{x}_{j}}\equiv \frac{{{I}_{j}}}{L}.\end{eqnarray} \tag{ 7 }$$

Of course, the particular system being studied might be such that the interval of support of x is bounded; in all cases however these bounds are state-independent and their presence does not affect the reasonings presented here.

It is important to remember that ${{\rho}_{\text{dist}}}$ carries the complete characteristics of a given eigenstate, and is therefore (implicitly) labeled by the set {I}. Note also that it is not a continuously differentiable function of x, but rather a distribution visualized as an unevenly-spaced but equal-height delta comb defined by the quantum numbers. It exactly (to all orders in 1/L) obeys the identities

$$\begin{eqnarray}&&{\int}^{}\text{d}x{{\rho}_{\text{dist}}}\left(x;{{\left\{I\right\}}_{N}}\right)=\frac{N}{L},\quad {\int}^{}\text{d}x{{\rho}_{\text{dist}}}\left(x;{{\left\{I\right\}}_{N}}\right)\,f(x)=\frac{1}{L}\underset{j=1}{\overset{N}{\sum}}f\left({{x}_{j}}\right)\end{eqnarray} \tag{ 8 }$$

in which the integral is taken over the whole support of $\rho (x)$ and f (x) is any continuous function.

The thermodynamic limit, denoted ${{\lim}_{\text{Th}}}$ in the following, is defined as the simultaneous limits $L\to \infty$, $N\to \infty$ with N/L fixed. For later purposes, we will also consider a more flexible viewpoint here and split the space of quantum numbers into 'boxes' (labeled by an integer i) of size l_i chosen such that ${{\lim}_{\text{Th}}}\frac{{{l}_{i}}}{L}=0$ and ${{\lim}_{\text{Th}}}{{l}_{i}}=+\infty$ $\forall i$, giving us a 'box regularized' thermodynamic limit which we will denote by ${{\lim}_{\text{Th},\text{reg}}}$.

For a given eigenstate, within box i, there will be a positive integer number n_i of occupied quantum numbers, allowing to define an average density ${{\rho}_{i}}$ in box i as

$$\begin{eqnarray}&&{{\rho}_{i}}=\frac{{{n}_{i}}}{{{l}_{i}}},\quad \text{with}\quad 0\leqslant {{\rho}_{i}}\leqslant 1~\forall i.\end{eqnarray} \tag{ 9 }$$

The specification of the box fillings ${{\rho}_{i}}$ does not label eigenstates uniquely. In order to do so, we also need to specify the particular configuration of quantum numbers within all boxes. Assuming that c_i (which is an index running over $\left(\begin{array}{c} {{l}_{i}} \\ {{n}_{i}} \end{array}\right)$ values) labels such a configuration in box i, we can thus uniquely label an eigenstate for a given box regularization by its set of density and in-box configuration labels, $|\left\{I\right\}\rangle \equiv |\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\rangle$, the resolution of the identity becoming

$$\begin{eqnarray}&&\mathbf{1}=\underset{\left\{I\right\}}{\sum}|\left\{I\right\}\rangle \langle \left\{I\right\}|=\underset{\left\{{{\rho}_{i}}\right\}}{\sum}\underset{\left\{{{c}_{i}}\right\}}{\sum}|\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\rangle \langle \left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}|.\end{eqnarray} \tag{ 10 }$$

In the Hilbert space, for a given box regularization, there will be many states with the same set of box densities $\left\{{{\rho}_{i}}\right\}$. This number is simply

$$\begin{eqnarray}&&\underset{\left\{{{c}_{i}}\right\}}{\sum}1{{|}_{\text{fixed}~\left\{{{\rho}_{i}}\right\}}}=\underset{i}{\prod}\left(\begin{array}{c} {{l}_{i}} \\ {{n}_{i}} \end{array}\right)\equiv {{\text{e}}^{S_{\left\{{{\rho}_{i}}\right\}}^{\text{YY}}}}\end{eqnarray} \tag{ 11 }$$

whose logarithm has a leading extensive term which is the well-known Yang–Yang entropy [25].

Within the set of states at fixed $\left\{{{\rho}_{i}}\right\}$, there exists one with a 'maximally flat' configuration where, in all boxes i, the n_i occupied quantum numbers are distributed uniformly over the l_i available spaces². We will label this as {c_i}  =  {0}. We can then write

$$\begin{eqnarray}&&{{\rho}_{\text{dist}}}\left(x;\left\{I\right\}\right)={{\rho}_{\text{dist}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\right)=\rho \left(x;\left\{{{\rho}_{i}}\right\}\right)+{{\rho}_{\text{ib}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\right)\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \rho \left(x;\left\{{{\rho}_{i}}\right\}\right)\equiv {{\rho}_{\text{dist}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{0\right\}\right), \\ {{\rho}_{\text{ib}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\right)\equiv {{\rho}_{\text{dist}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\right)-\rho \left(x;\left\{{{\rho}_{i}}\right\}\right) \end{array}\end{eqnarray} \tag{ 13 }$$

respectively represent the 'maximally flat' part of the distribution, and its in-box microscopic features. In the regularized thermodynamic limit, we can effectively make the replacements

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}\rho \left(x;\left\{{{\rho}_{i}}\right\}\right)={{\rho}_{\text{sm}}}(x)\in {{C}^{\infty}},\quad \underset{\text{Th},\text{reg}}{{\lim}}{{\rho}_{\text{ib}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\right)=0\end{eqnarray} \tag{ 14 }$$

as far as any in-box-configuration insensitive quantity is concerned, in the sense that given a continuously differentiable function f(x), we have

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\text{Th},\text{reg}}{{\lim}}{\int}^{}\text{d}x\rho \left(x;\left\{{{\rho}_{i}}\right\}\right)\,f(x)={\int}^{}\text{d}x{{\rho}_{\text{sm}}}(x)\,f(x), \\ \underset{\text{Th},\text{reg}}{{\lim}}{\int}^{}\text{d}x{{\rho}_{\text{ib}}}\left(x;\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\right)\,f(x)=0. \end{array}\end{eqnarray} \tag{ 15 }$$

The smooth, differentiable function ${{\rho}_{\text{sm}}}(x)$ will be called the root density in order to distinguish it from the neither smooth nor differentiable distribution ${{\rho}_{\text{dist}}}$. Another convenient function to introduce is the hole density function, which is the complement of ${{\rho}_{\text{dist}}}$ in the space of quantum numbers, and the total density ${{\rho}^{t}}=\rho +{{\rho}^{h}}$ whose smooth version coincides with unity,

$$\begin{eqnarray}&&\rho _{\text{sm}}^{t}(x)={{\rho}_{\text{sm}}}(x)+\rho _{\text{sm}}^{h}(x)=1.\end{eqnarray} \tag{ 16 }$$

The notations up to now have involved integrations over quantum number space. It is customary to also express integrals in rapidity space, extending the Bethe equations (5) to the functional equation

$$\begin{eqnarray}&&{{\theta}_{\text{kin}}}\left(\lambda \right)+{\int}^{}\text{d}{{\lambda}^{\prime}}~{{\theta}_{\text{scat}}}\left(\lambda -{{\lambda}^{\prime}}\right){{\rho}_{\text{sm}}}\left({{\lambda}^{\prime}}\right)=2\pi ~x\left(\lambda \right)\end{eqnarray} \tag{ 17 }$$

where the function $x\left(\lambda \right)$ is now viewed as labelling the state. This function then defines a mapping between the quantum number and rapidity spaces. The transformation rule of δ functions then allows to write the densities directly in rapidity space,

$$\begin{eqnarray}&&{{\rho}_{\text{sm}}}\left(\lambda \right)={{\rho}_{\text{sm}}}\left(x\left(\lambda \right)\right)\frac{\text{d}x\left(\lambda \right)}{\text{d}\lambda},\quad \rho _{\text{sm}}^{t}\left(\lambda \right)=\frac{\text{d}x\left(\lambda \right)}{\text{d}\lambda}.\end{eqnarray} \tag{ 18 }$$

A technical aspect to bear in mind is that the distributions so defined are not made of equal-height delta functions, but rather of state-dependent heights, which are implicit functions of the interaction-induced backflows when moving rapidities around. The quantum numbers are the truly statistically-independent parameters from which the trace over states should be defined.

For completeness, we recall a few basic notions of the algebraic Bethe Ansatz which are useful for the discussions here. Readers looking for a more systematic introduction to the subject should consult [21] and references therein.

At the center of the ABA is the notion of the R-matrix, which is a matrix defined in the tensor product of two auxiliary spaces ${{\mathcal{A}}_{1}}\otimes {{\mathcal{A}}_{2}}$, obeying the Yang–Baxter relation

$$\begin{eqnarray}&&{{R}_{12}}\left(\lambda,\mu \right){{R}_{13}}\left(\lambda,\nu \right){{R}_{23}}\left(\mu,\nu \right)={{R}_{23}}\left(\mu,\nu \right){{R}_{13}}\left(\lambda,\nu \right){{R}_{12}}\left(\lambda,\mu \right).\end{eqnarray} \tag{ 19 }$$

One then introduces a monodromy matrix $T\left(\lambda \right)$ living in the tensor product of auxiliary and Hilbert spaces and depending on a (in general complex-valued) spectral parameter λ, and obeying the intertwining relation

$$\begin{eqnarray}&&{{R}_{12}}\left(\lambda,\mu \right){{T}_{1}}\left(\lambda \right){{T}_{2}}\left(\mu \right)={{T}_{2}}\left(\mu \right){{T}_{1}}\left(\lambda \right){{R}_{12}}\left(\lambda,\mu \right).\end{eqnarray} \tag{ 20 }$$

In the simplest cases, the auxiliary space is two-dimensional and the monodromy matrix is represented as

$$\begin{eqnarray}T\left(\lambda \right)\equiv \left(\begin{array}{c} A\left(\lambda \right) & B\left(\lambda \right) \\ C\left(\lambda \right) & D\left(\lambda \right) \end{array}\right)\end{eqnarray} \tag{ 21 }$$

in which $A\left(\lambda \right),B\left(\lambda \right),C\left(\lambda \right),D\left(\lambda \right)$ are operators acting in Hilbert space. The transfer matrix τ, defined as the trace (in auxiliary space) of the monodromy matrix (the second equality here is special to the case of two-dimensional auxiliary space)

$$\begin{eqnarray}&&\tau \left(\lambda \right)=\text{T}{{\text{r}}_{\mathcal{A}}}T\left(\lambda \right)=A\left(\lambda \right)+D\left(\lambda \right)\end{eqnarray} \tag{ 22 }$$

then has the fundamental property of self-commutation at arbitrarily chosen spectral parameters,

$$\begin{eqnarray}&&\left[\tau \left(\lambda \right),\tau \left(\mu \right)\right]=0,~\forall ~\lambda,\mu ~\in ~\mathbb{C}.\end{eqnarray} \tag{ 23 }$$

This naturally allows to build a set of commuting quantum charges from so-called trace identities

$$\begin{eqnarray}&&{{Q}_{n}}\propto \frac{{{\text{d}}^{n}}}{\text{d}{{\lambda}^{n}}}\ln \tau \left(\lambda \right){{|}_{\lambda =\xi}}\end{eqnarray} \tag{ 24 }$$

in which the constant prefactors can be chosen at leisure, and ξ is an evaluation point chosen to make the charges local (the R-matrix then becoming of permutation form).

All fundamental commutation relations of the model are contained in the algebra of monodromy matrix operators A, B, C, D issuing from (20) once an R-matrix and a monodromy matrix have been defined. The algebra is quadratic in operators. For example, the simplest nontrivial structure for an R-matrix takes the shape

$$\begin{eqnarray}R\left(\lambda,\mu \right)=\left(\begin{array}{c} 1 & 0 & 0 & 0 \\ 0 & b\left(\lambda,\mu \right) & c\left(\lambda,\mu \right) & 0 \\ 0 & c\left(\lambda,\mu \right) & b\left(\lambda,\mu \right) & 0 \\ 0 & 0 & 0 & 1 \end{array}\right)\end{eqnarray} \tag{ 25 }$$

with some constraints (which we do not write down here) on the b, c functions to ensure algebraic consistency through the Yang–Baxter relation. The intertwining relations (20) then imply

$$\begin{eqnarray}&&\left[A\left(\lambda \right),A\left(\mu \right)\right]=0,\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\left[B\left(\lambda \right),B\left(\mu \right)\right]=0,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&A\left(\lambda \right)B\left(\mu \right)=\frac{1}{b\left(\mu,\lambda \right)}B\left(\mu \right)A\left(\lambda \right)-\frac{c\left(\mu,\lambda \right)}{b\left(\mu,\lambda \right)}B\left(\lambda \right)A\left(\mu \right),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\left[B\left(\lambda \right),C\left(\mu \right)\right]=\frac{c\left(\lambda,\mu \right)}{b\left(\lambda,\mu \right)}\left(D\left(\mu \right)A\left(\lambda \right)-D\left(\lambda \right)A\left(\mu \right)\right),\end{eqnarray} \tag{ 29 }$$

and 12 other similar relations which we leave out here for the sake of brevity (see e.g. [21] for additional details).

Assuming the existence of a reference state $|0\rangle$ such that

$$\begin{eqnarray}&&A\left(\lambda \right)|0\rangle =a\left(\lambda \right)|0\rangle,\quad D\left(\lambda \right)|0\rangle =d\left(\lambda \right)|0\rangle,\quad C\left(\lambda \right)|0\rangle =0\end{eqnarray} \tag{ 30 }$$

where the functions $a\left(\lambda \right)$ and $d\left(\lambda \right)$ can be chosen arbitrarily (they should be viewed as model-defining choices), states in Hilbert space can be constructed by acting with the 'raising' operators $B\left(\lambda \right)$ according to

$$\begin{eqnarray}&&|{{\left\{{{\lambda}_{j}}\right\}}_{M}}\rangle \equiv \underset{j=1}{\overset{M}{\prod}}B\left({{\lambda}_{j}}\right)|0\rangle \end{eqnarray} \tag{ 31 }$$

for generic M and ${{\left\{{{\lambda}_{j}}\right\}}_{M}}$. Note the very important fact that the order in the product is immaterial, in view of the commutation relation (27). So written, the states are also not normalized. Again using the monodromy matrix element commutation relations, one can show that these states are eigenstates of the transfer matrix (22) provided the set of rapidities $\left\{\lambda \right\}$ obeys the Bethe equations

$$\begin{eqnarray}&&\frac{a\left({{\lambda}_{j}}\right)}{d\left({{\lambda}_{j}}\right)}\underset{l\ne j}{\prod}\frac{b\left({{\lambda}_{j}},{{\lambda}_{l}}\right)}{b\left({{\lambda}_{l}},{{\lambda}_{j}}\right)}=1~ ~\forall ~j=1,...,M,\end{eqnarray} \tag{ 32 }$$

whose logarithm (allowing the introduction of quantum numbers {I}, which allow to classify eigenstates once the Pauli principle of non-coincident rapidities has been implemented) we have written as (5). These states diagonalize the transfer matrix

$$\begin{eqnarray}&&\tau \left(\lambda \right)|{{\left\{{{\lambda}_{j}}\right\}}_{M}}\rangle =\tau \left(\lambda |{{\left\{{{\lambda}_{j}}\right\}}_{M}}\right)|{{\left\{{{\lambda}_{j}}\right\}}_{M}}\rangle,\end{eqnarray} \tag{ 33 }$$

with eigenvalues

$$\begin{eqnarray}&&\tau \left(\lambda |{{\left\{{{\lambda}_{j}}\right\}}_{M}}\right)=a\left(\lambda \right)\underset{j=1}{\overset{M}{\prod}}{{b}^{-1}}\left({{\lambda}_{j}},\lambda \right)+d\left(\lambda \right)\underset{j=1}{\overset{M}{\prod}}{{b}^{-1}}\left(\lambda,{{\lambda}_{j}}\right),\end{eqnarray} \tag{ 34 }$$

this defining the eigenvalue of all charges (24).

A further pillar of the algebraic Bethe Ansatz is Slavnov's theorem [26], which gives the overlap between states

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}_{M}}\left(\left\{\mu \right\},\left\{\lambda \right\}\right)=\langle 0|\underset{j=1}{\overset{M}{\prod}}C\left({{\mu}_{j}}\right)\underset{k=1}{\overset{M}{\prod}}B\left({{\lambda}_{k}}\right)|0\rangle \end{array}\end{eqnarray} \tag{ 35 }$$

as a computationally convenient determinant, provided either $\left\{\lambda \right\}$ or $\left\{\mu \right\}$ obeys the Bethe equations:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}_{M}}\left(\left\{\mu \right\},\left\{\lambda \right\}\right)=\frac{\prod\nolimits_{j=1}^{M}\prod\nolimits_{k=1}^{M}\varphi \left({{\mu}_{j}}-{{\lambda}_{k}}\right)}{\prod\nolimits_{j<k}\varphi \left({{\mu}_{j}}-{{\mu}_{k}}\right)\prod\nolimits_{j>k}\varphi \left({{\lambda}_{j}}-{{\lambda}_{k}}\right)}\det T\left(\left\{\mu \right\},\left\{\lambda \right\}\right), \end{array}\end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray}&&{{T}_{ab}}=\frac{\partial}{\partial {{\lambda}_{a}}}\tau \left({{\mu}_{b}}|\left\{\lambda \right\}\right)\end{eqnarray} \tag{ 37 }$$

and φ is a model-dependent scalar function. The importance of Slavnov's theorem cannot be overemphasized, since it allows for the computation of matrix elements in integrable models once the solution to the 'quantum inverse problem' (the mapping of physical operators to monodromy matrix operators, [27–33]) is known.

The norm of an eigenstate can similarly be computed algebraically using the commutation relations between monodromy matrix entries, or as the appropriate limit of Slavnov's theorem, yielding the celebrated Gaudin formula [20, 34, 35].

Let us go back to our original problem and consider an arbitrary wavefunction at t  =  0. We assume that from this instant onwards, the time evolution occurs according to the Bethe Ansatz-solvable Hamiltonian H whose normalized eigenstates are labeled by sets of quantum numbers {I }. The initial state is exactly decomposed in the basis of these eigenstates according to

$$\begin{eqnarray}&&|\Psi(t=0)\rangle =\underset{\left\{I\right\}}{\sum}{{\text{e}}^{-S_{\left\{I\right\}}^{\Psi}}}|\left\{I\right\}\rangle \end{eqnarray} \tag{ 38 }$$

where we have defined the overlap coefficients

$$\begin{eqnarray}&&S_{\left\{I\right\}}^{\Psi}=-\ln \langle \left\{I\right\}|\Psi(t=0)\rangle ~\in ~\mathbb{C}.\end{eqnarray} \tag{ 39 }$$

One crucial property is that since we are working with normalized states, the real parts of the overlap coefficients are bounded from below,

$$\begin{eqnarray}&&\exists ~{{S}_{\text{min}}}\in \mathbb{R}\geqslant 0~|~\Re S_{\left\{I\right\}}^{\Psi}\geqslant {{S}_{\text{min}}}~\forall \left\{I\right\}\end{eqnarray} \tag{ 40 }$$

and tend to infinity for states with vanishing overlap.

In the eigenbasis, the time-dependent Schrödinger equation is now trivially solved by using $H|\left\{I\right\}\rangle ={{\omega}_{\left\{I\right\}}}|\left\{I\right\}\rangle$, allowing us to write the exact time-dependent wavefunction as

$$\begin{eqnarray}&&|\Psi(t)\rangle =\underset{\left\{I\right\}}{\sum}{{\text{e}}^{-S_{\left\{I\right\}}^{\Psi}-\text{i}{{\omega}_{\left\{I\right\}}}t}}|\left\{I\right\}\rangle.\end{eqnarray} \tag{ 41 }$$

The expectation value (2) which we are centrally interested in can then of course be written, without any approximation, as a double Hilbert space-sized summation

$$\begin{eqnarray}&&\bar{\mathcal{O}}(t)=\frac{\sum\nolimits_{\{{{I}^{l}}\}}\sum\nolimits_{\{{{I}^{r}}\}}{{\text{e}}^{-{{(S_{\{{{I}^{l}}\}}^{\Psi})}^{\ast}}-S_{\{{{I}^{r}}\}}^{\Psi}+\text{i}({{\omega}_{\{{{I}^{l}}\}}}-{{\omega}_{\{{{I}^{r}}\}}})t}}\langle \{{{I}^{l}}\}|\mathcal{O}|\{{{I}^{r}}\}\rangle}{\sum\nolimits_{\{I\}}{{\text{e}}^{-2\Re S_{\{I\}}^{\Psi}}}}.\end{eqnarray} \tag{ 42 }$$

Such a double summation is in general too difficult to perform due to the exceedingly large number of terms it contains. Our aim here is to show that in the thermodynamic limit, under mild assumptions, this summation can be drastically simplified.

Let us begin by looking at the denominator of (2) or (42), namely the (time-independent) normalization of the initial state:

$$\begin{eqnarray}&&\langle \Psi(t)|\Psi(t)\rangle =\underset{\left\{I\right\}}{\sum}{{\text{e}}^{-2\Re S_{\left\{I\right\}}^{\Psi}}}.\end{eqnarray} \tag{ 43 }$$

Using our resolution of the identity (10), we can write this (using a self-explanatory notation) as

$$\begin{eqnarray}&&\langle \Psi(t)|\Psi(t)\rangle =\underset{\left\{{{\rho}_{i}}\right\}}{\sum}\underset{\left\{{{c}_{i}}\right\}}{\sum}{{\text{e}}^{-2\Re S_{\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}}^{\Psi}}}.\end{eqnarray} \tag{ 44 }$$

If we are using normalized states at all steps, this summation is of course equal to one. The detailed distribution of contributions is however not trivial, and the weight distribution is utterly undemocratic. For a generic initial state, the overlaps can be rapidly-varying numbers as we move around the Hilbert space. Modifying the quantum number of a single particle even by the smallest allowable unit can for example change the overlap by a factor of order one. In fact, the overlaps can even suddenly identically vanish if some discrete symmetry requirement is violated. Assuming all such discrete symmetries have been handled by a proper partitioning of the Hilbert space³, the remaining overlaps do not by any means have to be smooth functions of the state's quantum numbers. On the other hand, one can expect that the extensive part of the logarithm of the overlaps is insensitive to the microscopic details of the state, and rather depends only on the overall root distribution. To express this formally, for a given set of box fillings $\left\{{{\rho}_{i}}\right\}$, we can define an effective box-averaged overlap according to

$$\begin{eqnarray}&&{{\text{e}}^{-S_{\left\{{{\rho}_{i}}\right\}}^{o}}}\equiv {{\text{e}}^{-S_{\left\{{{\rho}_{i}}\right\}}^{\text{YY}}}}\underset{\left\{{{c}_{i}}\right\}}{\sum}{{\text{e}}^{-2\Re S_{\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}}^{\Psi}}}\end{eqnarray} \tag{ 45 }$$

in which $S_{\left\{{{\rho}_{i}}\right\}}^{o}$ is a real-valued function of the fillings $\left\{{{\rho}_{i}}\right\}$, representing the overlaps, which becomes a well-defined, differentiable functional of the smooth distribution ${{\rho}_{\text{sm}}}$ (14) in the thermodynamic limit,

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}S_{\left\{{{\rho}_{i}}\right\}}^{o}\equiv {{S}^{o}}\,\left[{{\rho}_{\text{sm}}}\right].\end{eqnarray} \tag{ 46 }$$

The same smoothness in the thermodynamic limit is of course automatically true (and traditionally tacitly assumed) of the Yang–Yang entropy,

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}S_{\left\{{{\rho}_{i}}\right\}}^{\text{YY}}\equiv {{S}^{\text{YY}}}\left[{{\rho}_{\text{sm}}}\right].\end{eqnarray} \tag{ 47 }$$

Once the in-box summations have been performed in this way, the remaining summation over box fillings can be interpreted as a conventional functional integral over continuously differentiable functions ${{\rho}_{\text{sm}}}$,

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}\underset{\left\{{{\rho}_{i}}\right\}}{\sum}\left(...\right)={{{\int}^{}}_{{{\rho}_{\text{sm}}}\in {{C}^{\infty}}}}D{{\rho}_{\text{sm}}}\left(...\right)\end{eqnarray} \tag{ 48 }$$

(since we are still working with distributions in quantum number space, where statistical independence holds, there is no extra Jacobian in this functional integral) so the normalization summation (44) can be rewritten as

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}\langle \Psi(t)|\Psi(t)\rangle ={\int}^{}D{{\rho}_{\text{sm}}}~{{\text{e}}^{-{{S}^{Q}}\left[{{\rho}_{\text{sm}}}\right]}}\end{eqnarray} \tag{ 49 }$$

where the 'Quench Action' functional is defined as the difference between the pseudo-energy obtained from the overlaps and the Yang–Yang entropy of the state,

$$\begin{eqnarray}&&{{S}^{Q}}\left[\rho \right]={{S}^{o}}\left[\rho \right]-{{S}^{\text{YY}}}\left[\rho \right].\end{eqnarray} \tag{ 50 }$$

By construction, the Quench Action thus represents a sort of equivalent of a free energy for out-of-equilibrium situations. It is an extensive (any non-extensive part takes on large positive values, and thus does not influence the expectation values considered), real-valued functional which is bounded from below due to the state normalization constraint. One of the appeals of the Quench Action approach is that it is amenable to all the standard field-theoretical methods we are accustomed to in (quantum) statistical mechanics. This will be discussed in more detail later on in the perspectives.

As a first step towards a more useful representation, let us apply a steepest-descent reasoning in the evaluation of (49). Assuming that there exists a single minimum⁴, the functional integral (49) can be evaluated in a saddle-point approximation (the applicability of the saddle-point logic resting on the system size going to infinity),

$$\begin{eqnarray}&&{\int}^{}D{{\rho}_{\text{sm}}}~{{\text{e}}^{-{{S}^{Q}}\left[{{\rho}_{\text{sm}}}\right]}}={{\text{e}}^{-{{S}^{Q}}\left[{{\rho}_{\text{sp}}}\right]}}~\text{Det}{{\left[\frac{T}{2\pi}\right]}^{-1/2}}\left(1+{{C}_{L}}\right)\end{eqnarray} \tag{ 51 }$$

where C_L represent corrections vanishing in the thermodynamic limit, and $T\left(\lambda,\mu \right)=$ $\left. \frac{{{\delta}^{2}}{{S}^{Q}}\left[\rho \right]}{\delta \rho \left(\lambda \right)\delta \rho \left(\mu \right)}\right|_{{{\rho}_{\text{sp}}}}$ is the functional Hessian of the Quench Action evaluated at the (continuously differentiable) saddle-point distribution ${{\rho}_{\text{sp}}}$, which is determined as the distribution that satisfies the generalized TBA (GTBA, see for example [37, 38]) equilibrium condition

$$\begin{eqnarray}&&\left. \frac{\delta {{S}^{Q}}\left[\rho \right]}{\delta \rho}\right|_{{{\rho}_{\text{sp}}}}=0.\end{eqnarray} \tag{ 52 }$$

The Hessian (provided it doesn't vanish) gives only subleading contributions in the thermodynamic limit, and can be neglected if one focuses on the dominant contributions only. Including subdominant terms is procedurally straightforward from the construction above.

In practice, the equations derived from the saddle-point condition (52) are morphologically identical to the thermodynamic Bethe Ansatz equations one obtains when treating finite-temperature equilibrium integrable models. In TBA (see [25, 22] and references therein), one indeed performs a saddle-point analysis, which becomes exact in the thermodynamic limit, the functional integral weight being simply the (exponential of minus the) free energy of the system. The Quench Action (50) however contains terms which depend on the density function ρ in ways distinct from the thermal free energy, meaning that the 'driving terms' in equations (52) are distinct from the usual TBA ones. Explicit examples of these will be given in the next section.

Now is a good time to discuss a subtle point: is the Quench Action (50) equivalent to the GGE? The Quench Action is by construction a mathematically exact representation of the diagonal ensemble; the GGE, when implemented correctly, converges to it in the thermodynamic limit. If the saddle-point approximation can be used in both cases to find a steady state, then these steady states must be the same. In this sense, equivalence is obvious. That said, the fact that the QA and the GGE share the saddle point does not mean that beyond-saddle-point features must coincide. More importantly, the point is that the QA and the GGE are founded on different footholds, meaning that approximations within one approach will not be expressible within the language of the other. For example, approximating the GGE by truncating the set of charges used has no meaningful equivalent within the QA approach. Conversely, a clever approximation scheme for evaluating the extensive parts of the overlap logarithms has no obvious translation to the language of conserved charges. Thus, though properly-implemented and performed calculations within the QA and GGE schemes should agree on the steady state, the two approaches remain quite distinct on the practical level. A more fundamental difference between the approaches is however discussed in the next section.

Let us now tackle the more challenging (but physically rich and interesting) question of the time evolution. In the context of integrable models, much work has already been performed on this issue, see for example [39–58].

Within the Quench Action formalism, obtaining the time dependence involves evaluation of the numerator of (42),

$$\begin{eqnarray}&&\langle \Psi(t)|\mathcal{O}|\Psi(t)\rangle =\underset{\{{{I}^{l}}\}}{\sum}\underset{\{{{I}^{r}}\}}{\sum}{{\text{e}}^{-{{(S_{\{{{I}^{l}}\}}^{\Psi})}^{\ast}}-S_{\{{{I}^{r}}\}}^{\Psi}}}{{\text{e}}^{\text{i}({{\omega}_{\{{{I}^{l}}\}}}-{{\omega}_{\{{{I}^{r}}\}}})t}}\langle \{{{I}^{l}}\}|\mathcal{O}|\{{{I}^{r}}\}\rangle.\end{eqnarray} \tag{ 53 }$$

In the thermodynamic limit, we expect that only a minority of states significantly contribute to this double sum. A major simplification comes from exploiting the fact that for a typical physical operator $\mathcal{O}$ (local density, particle addition/removal, etc), the matrix elements $\langle \{{{I}^{l}}\}|\mathcal{O}|\{{{I}^{r}}\}\rangle$ are rapidly-decreasing functions of the 'difference' between the bra and ket states, in other words of the number of displaced quantum numbers from one state to the other. We will call the operator $\mathcal{O}$ a 'weak' operator if its matrix elements $\langle \{{{I}^{l}}\}|\mathcal{O}|\{{{I}^{r}}\}\rangle$ are negligible unless $\{{{I}^{l}}\}=\{{{I}^{r}}\}+\text{excitations}~\text{carrying}~\text{zero}~\text{entropy}.$ In other words, a weak operator $\mathcal{O}$ does not produce finite-entropy modifications of the state it acts upon. An equivalent way of saying this is that the linear decomposition

$$\begin{eqnarray}&&\mathcal{O}|\{I\}\rangle =\underset{\mathbf{e}}{\sum}{{o}_{\{I\};\mathbf{e}}}|\{I\};\mathbf{e}\rangle,\quad {{o}_{\{I\};\mathbf{e}}}\in \mathbb{C}\end{eqnarray} \tag{ 54 }$$

representing the action of $\mathcal{O}$ on a certain state can be truncated to a discrete, sub-entropically large number of excitations as one approaches the thermodynamic limit, while obeying all available sum rules to arbitrary accuracy. To put it quantitatively, let us use the operator $\mathcal{O}$'s normalized elements

$$\begin{eqnarray}&&{{n}_{\left\{I\right\};\mathbf{e}}}\equiv {{o}_{\left\{I\right\};\mathbf{e}}}{{\left[\underset{\mathbf{e}}{\sum}|{{o}_{\left\{I\right\};\mathbf{e}}}{{|}^{2}}\right]}^{-1/2}}\end{eqnarray} \tag{ 55 }$$

to define a Shannon-like operator entropy

$$\begin{eqnarray}&&S_{\{I\}}^{\mathcal{O}}\equiv -\underset{\mathbf{e}}{\sum}(|{{n}_{\{I\};\mathbf{e}}}{{|}^{2}}\ln |{{n}_{\{I\};\mathbf{e}}}{{|}^{2}}+\,(1\,-\,|{{n}_{\{I\};\mathbf{e}}}{{|}^{2}})\ln (1\,-\,|{{n}_{\{I\};\mathbf{e}}}{{|}^{2}})).\end{eqnarray} \tag{ 56 }$$

A weak operator is then an operator which has a subextensive operator entropy. Note that this statement depends on both the operator and on the state on which it is applied, and is moreover basis-dependent: the operator entropy is calculated in the eigenstates basis (since we are interested in dephasings under time evolution). It implicitly depends on the correspondence between the structure of the constituents of $\mathcal{O}$ (for example, $\mathcal{O}$ could be a simple sum of local operators) and the operators corresponding to the creation of eigenstates (namely, products of B operators). If these are similar, the operator entropy is low. One also generally expects a finite product of weak operators to itself be weak. Note that multiple-point, time-split operator insertions like ${{\text{e}}^{\text{i}({{t}_{a}}-{{t}_{b}})H}}{{\mathcal{O}}_{a}}{{\text{e}}^{-\text{i}({{t}_{a}}-{{t}_{b}})H}}{{\mathcal{O}}_{b}}$ representing dynamical correlations, are also weak operators if ${{\mathcal{O}}_{a}}$ and ${{\mathcal{O}}_{b}}$ are weak.

Assuming that $\mathcal{O}$ is a weak operator, the numerator of (42) simplifies to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \Psi(t)|\mathcal{O}|\Psi(t)\rangle =\underset{\{{{\rho}_{i}}\}}{\sum}\underset{\{{{c}_{i}}\}}{\sum}\underset{\mathbf{e}}{\sum}{{\text{e}}^{-{{(S_{\{{{\rho}_{i}}\};\{{{c}_{i}}\}}^{\Psi})}^{\ast}}-S_{\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}}^{\Psi}}}\times {{\text{e}}^{\text{i}({{\omega}_{\{{{\rho}_{i}}\};\{{{c}_{i}}\}}}-{{\omega}_{\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}}})t}}\langle \{{{\rho}_{i}}\};\{{{c}_{i}}\}|\mathcal{O}|\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}\rangle. \end{array}\end{eqnarray} \tag{ 57 }$$

To proceed further, we note that the thermodynamically finite energy difference between the bra and ket states depends only on the discrete excitations and the smooth density ${{\rho}_{\text{sm}}}$, but not on the in-box configurations {c_i}. This is a very general and robust feature of states in the thermodynamic limit. We will write

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\text{Th},\text{reg}}{{\lim}}({{\omega}_{\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}}}-{{\omega}_{\{{{\rho}_{i}}\};\{{{c}_{i}}\}}})\equiv {{\omega}_{\mathbf{e}}}\,[{{\rho}_{\text{sm}}}]. \end{array}\end{eqnarray} \tag{ 58 }$$

Since the energy- and time-dependent phase can now be factorized out of the {c_i} configuration sums, we can define box-averaged combinations of the overlaps and matrix element products. To do this, we assume that the weak operator $\mathcal{O}$ we consider is such that its matrix elements are invariant under a simultaneous shift of in-box configurations in the bra and ket states, up to corrections which vanish in the thermodynamic limit. We will then call the operator 'smooth'. Symbolically, we can write that a smooth operator $\mathcal{O}$ is such that

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}\frac{\langle \{{{\rho}_{i}}\};\{c_{i}^{\prime}\}|\mathcal{O}|\{{{\rho}_{i}}\};\{c_{i}^{\prime}\};\mathbf{e}\rangle}{\langle \{{{\rho}_{i}}\};\{{{c}_{i}}\}|\mathcal{O}|\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}\rangle}=1+o(1)\quad \forall ~\{{{c}_{i}}\},\{c_{i}^{\prime}\}.\end{eqnarray} \tag{ 59 }$$

Equivalently, a smooth operator $\mathcal{O}$ commutes with operators enforcing only in-box modifications of the quantum numbers, up to terms vanishing in the regularized thermodynamic limit. If we define in-box shuffling operators

$$\begin{eqnarray}&&{{P}_{\left\{{{c}_{i}}\right\}\left\{c_{i}^{\prime}\right\}}}\equiv |\left\{{{c}_{i}}\right\}\rangle \langle \left\{c_{i}^{\prime}\right\}|,\end{eqnarray} \tag{ 60 }$$

the statement that operator $\mathcal{O}$ is smooth is equivalent to the vanishing of the commutation of $\mathcal{O}$ with all in-box shufflings, as far as all expectation values are concerned,

$$\begin{eqnarray}&&\langle \{{{\rho}_{i}}\}|[\mathcal{O},{{P}_{\{{{c}_{i}}\}\{c_{i}^{\prime}\}}}]|\{{{\rho}_{i}}\}\rangle =o(1)\quad \forall ~\{{{c}_{i}}\},\{c_{i}^{\prime}\}.\end{eqnarray} \tag{ 61 }$$

The assumption that $\mathcal{O}$ is a weak, smooth operator allows us to considerably simplify our equations by picking a particular in-box configuration as a representative and to take the matrix element out of the in-box configuration sum, writing

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\{{{c}_{i}}\}}{\sum}{{\text{e}}^{-{{(S_{\{{{\rho}_{i}}\};\{{{c}_{i}}\}}^{\Psi})}^{\ast}}-S_{\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}}^{\Psi}}}\langle \{{{\rho}_{i}}\};\{{{c}_{i}}\}|\mathcal{O}|\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}\rangle \equiv {{\text{e}}^{-S_{\{{{\rho}_{i}}\}}^{o}+S_{\{{{\rho}_{i}}\}}^{\text{YY}}}}\,{{\text{e}}^{-\delta {{S}_{\{{{\rho}_{i}}\};\mathbf{e}}}}}\langle \{{{\rho}_{i}}\};\{0\}|\mathcal{O}|\{{{\rho}_{i}}\};\{0\};\mathbf{e}\rangle \end{array}\end{eqnarray} \tag{ 62 }$$

where the box average of the overlap products

$$\begin{eqnarray}&&{{\text{e}}^{-\delta {{S}_{\{{{\rho}_{i}}\};\mathbf{e}}}}}\equiv {{\text{e}}^{S_{\{{{\rho}_{i}}\}}^{o}-S_{\{{{\rho}_{i}}\}}^{\text{YY}}}}\underset{\{{{c}_{i}}\}}{\sum}{{\text{e}}^{-{{(S_{\{{{\rho}_{i}}\};\{{{c}_{i}}\}}^{\Psi})}^{\ast}}-S_{\{{{\rho}_{i}}\};\{{{c}_{i}}\};\mathbf{e}}^{\Psi}}}\end{eqnarray} \tag{ 63 }$$

is well-defined and naturally becomes a smooth functional of ${{\rho}_{\text{sm}}}$ and function of the discrete excitations in the regularized thermodynamic limit,

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}\delta {{S}_{\{{{\rho}_{i}}\};\mathbf{e}}}\equiv \delta {{S}_{\mathbf{e}}}[{{\rho}_{\text{sm}}}]\end{eqnarray} \tag{ 64 }$$

with condition $\delta {{S}_{0}}[\rho ]=0$ by definition. Using this 'differential overlap' $\delta S$, the thermodynamic limit of our expression for the time-dependent correlator thus becomes, up to vanishing corrections,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\text{Th},\text{reg}}{{\lim}}\bar{\mathcal{O}}(t)=\frac{{\int}^{}\mathcal{D}{{\rho}_{\text{sm}}}{{\text{e}}^{-{{S}^{Q}}[{{\rho}_{\text{sm}}}]}}~\lim\nolimits_{\text{Th},\text{reg}}\sum\nolimits_{\mathbf{e}}{{\text{e}}^{-\delta {{S}_{\mathbf{e}}}[{{\rho}_{\text{sm}}}]-\text{i}{{\omega}_{\mathbf{e}}}[{{\rho}_{\text{sm}}}]t}}\langle {{\rho}_{\text{sm}}}|\mathcal{O}|{{\rho}_{\text{sm}}};\mathbf{e}\rangle}{{\int}^{}\mathcal{D}{{\rho}_{\text{sm}}}~{{\text{e}}^{-{{S}^{Q}}[{{\rho}_{\text{sm}}}]}}} \end{array}\end{eqnarray} \tag{ 65 }$$

where we interpret the summation in the numerator as the thermodynamic limit of a summation regularized by performing it at a large but not infinite size (in order to keep the matrix elements individually finite). $\langle {{\rho}_{\text{sm}}}|\mathcal{O}|{{\rho}_{\text{sm}}};\mathbf{e}\rangle$ is thus viewed at this stage as the matrix element obtained in a arbitrary but fixed regularization of ${{\rho}_{\text{sm}}}$ in terms of a box-regularized $\{{{\rho}_{i}}\}$ and a fixed choice of in-box configurations {c_i}, for example $\langle \{{{\rho}_{i}}\};\{0\}|\mathcal{O}|\{{{\rho}_{i}}\};\{0\};\mathbf{e}\rangle$. These matrix elements are completely defined (exactly, irrespective of the interaction parameter, to all orders in inverse system size) by their algebraic Bethe Ansatz representation [28, 29, 59–68] (typically in terms of a determinant or summations thereof), and thus completely free of singularities. They are however very strongly dependent on the relative microscopic positioning of the quantum numbers between the bra and ket states, and thus cannot be viewed as a 'smooth', slowly varying function of the relevant parameters here. The summation over $\mathbf{e}$ should be interpreted at this stage as running over all excitations defined by allowable changes in the quantum number configuration for a given microscopic regularization, which essentially splits up into in-box entropy-like summations and dispersing excitation summation. This last summation itself has however a completely well-defined thermodynamic limit. We discuss this point in more detail later on in section 3.4.

Formula (65) is already considerably simplified as compared to the original starting point (2), in the sense that only a single functional integral remains. It is of quite general applicability, being valid for arbitrary times t and any smooth operator $\mathcal{O}$. Under mild assumptions, we can however proceed further using a saddle-point evaluation. We can assume that the differential overlap $\delta {{S}_{\mathbf{e}}}\left[{{\rho}_{\text{sm}}}\right]$, encoding the correlation between quench overlaps of nearby states in Hilbert space, is non-extensive and thus does not (to leading order in system size) shift the saddle-point of the numerator as compared to that of the denominator. Note that this assumption is completely natural: modifying the quantum numbers of a single particle typically modifies the overlap by a factor which is algebraic rather than exponential in system size. Then, for operators $\mathcal{O}$ whose matrix elements are not exponentially large in system size (and thus also do not shift the saddle point; we will then call $\mathcal{O}$ thermodynamically finite), we can use the previously-obtained saddle point ${{\rho}_{\text{sp}}}$ (52) and get the yet further significantly simplified form (explicitly reintroducing for convenience the symmetry between putting the excitations in the bra or the ket, which is implicitly present but not manifest in (65))

$$\begin{eqnarray}\begin{array}{*{35}{l}} \underset{\text{Th},\text{reg}}{{\lim}}\bar{\mathcal{O}}(t) & =\underset{\text{Th},\text{reg}}{{\lim}}\frac{1}{2}\underset{\mathbf{e}}{\sum}[{{\text{e}}^{-\delta {{S}_{\mathbf{e}}}[{{\rho}_{\text{sp}}}]-\text{i}{{\omega}_{\mathbf{e}}}[{{\rho}_{\text{sp}}}]t}}\langle {{\rho}_{\text{sp}}}|\mathcal{O}|{{\rho}_{\text{sp}}};\mathbf{e}\rangle +\,{{\text{e}}^{-\delta S_{\mathbf{e}}^{\ast}[{{\rho}_{\text{sp}}}]+\text{i}{{\omega}_{\mathbf{e}}}[{{\rho}_{\text{sp}}}]t}}\langle {{\rho}_{\text{sp}}};\mathbf{e}|\mathcal{O}|{{\rho}_{\text{sp}}}\rangle ]. \end{array}\end{eqnarray} \tag{ 66 }$$

This equation is perhaps the most crucial formula of the Quench Action formalism. It is applicable provided operator $\mathcal{O}$ is

• (i)  
  weak,
• (ii)  
  smooth,
• (iii)  
  thermodynamically finite

(in particular, note that these conditions are fulfilled quite generally by local operators). Note the crucial fact that we have not assumed large times while deriving this formula; in fact, we conjecture that expression (66), evaluated in the thermodynamic limit, is valid for all times t  >  0 provided the operator $\mathcal{O}$ is weak, smooth and thermodynamically finite. At large times, meaning at times much larger than the lowest considered excitation energy (that being of order 1/L, and not ${{\text{e}}^{-(cst)L}}$ like the mean energy spacing), the formula simplifies to its diagonal part

$$\begin{eqnarray}&&\underset{\text{Th},\text{reg}}{{\lim}}\bar{\mathcal{O}}(t)=\underset{\text{Th},\text{reg}}{{\lim}}\langle {{\rho}_{\text{sp}}}|\mathcal{O}|{{\rho}_{\text{sp}}}\rangle \end{eqnarray} \tag{ 67 }$$

with right-hand side evaluated on any microscopic realization of the saddle-point state. This parallels the microcanonical sum used in [69], which can be interpreted as a generalization of the Eigenstate Thermalization Hypothesis [18, 70–74].

In terms of computational complexity, the representation (66) thus displays the fact that the Quench Action 'starts from the steady state' and works its way backwards in time upon the addition of more and more excitations around the steady state.

Yet a further slightly simplified version of (66) can be obtained by exploiting the fact that energies of the individual discrete excitations around the saddle point are additive. This is a well-known general feature of the thermodynamic limit of the exact eigenstates we are playing with. More precisely, if $\mathbf{e}$ represents the set of excitations (particles, holes, etc) which have been created on the state described by the thermodynamic root distribution ${{\rho}_{\text{sp}}}$, we have

$$\begin{eqnarray}&&{{\omega}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]=\underset{i}{\sum}\varepsilon \left({{e}_{i}}\right)\end{eqnarray} \tag{ 68 }$$

in which $\varepsilon \left({{e}_{i}}\right)$ is a single function representing the dispersion relation for excitations around the saddle point (and obtainable via the GTBA).

We can even go one step further and assume that, similarly to the energies, the differential overlap splits into decoupled functions of the individual excitations⁵, namely that there exists a single complex-valued function s such that

$$\begin{eqnarray}&&\delta {{S}_{\mathbf{e}}}\left[\rho \right]=\underset{i}{\sum}s\left({{e}_{i}}\right)+...\end{eqnarray} \tag{ 69 }$$

encoding the overlap differences between states appearing in the sum (66). Such a function can be read from considering the scaling of the exact overlaps as system size goes to infinity, an example being formula (109) (from (A14) of [36]). In fact, one cannot a priori exclude the possibility that (69) can be oversimplified, and that one should include a whole series of additional many-body terms

$$\begin{eqnarray}&&\delta {{S}_{\mathbf{e}}}\left[\rho \right]=\underset{i}{\sum}s\left({{e}_{i}}\right)+\underset{{{i}_{1}}<{{i}_{2}}}{\sum}{{s}^{(2)}}\left({{e}_{{{i}_{1}}}},{{e}_{{{i}_{2}}}}\right)+...\end{eqnarray} \tag{ 70 }$$

In most situations however the two-body term (and higher ones) would carry only finite-size corrections. The fact remains that the function s (and eventually it's many-body corrections) can be obtained directly from the exact overlaps. We will provide examples of this in the next section.

We can in fact push the preceding logic one step further. Let us examine in some more detail the required matrix elements

$$\begin{eqnarray}&&\langle {{\rho}_{\text{sp}}}|\mathcal{O}|{{\rho}_{\text{sp}}};\mathbf{e}\rangle.\end{eqnarray} \tag{ 71 }$$

Physically relevant operators are such that the value of such a matrix element decreases upon inserting more and more excitations. Moreover, the factor by which the matrix element decreases upon adding an excitation, depends primarily on the added excitation, with additional corrections coming from the interplay between excitations. This leads us to propose to represent the matrix elements (71) in terms of weights h as

$$\begin{eqnarray}&&\langle {{\rho}_{\text{sp}}}|\mathcal{O}|{{\rho}_{\text{sp}}};\mathbf{e}\rangle ={{\text{e}}^{-h_{{{\rho}_{\text{sp}}};\mathbf{e}}^{\mathcal{O}}}},\end{eqnarray} \tag{ 72 }$$

these weights taking complex values in general; for a normalized operator, the real part is however naturally bounded from below. One can of course 'lift' this representations to the operator level, writing $\mathcal{O}={{\text{e}}^{-\hat{h}_{{{\rho}_{\text{sp}}}}^{\mathcal{O}}}}$ in terms of the 'operator pseudo-Hamiltonian' $\hat{h}_{{{\rho}_{\text{sp}}}}^{\mathcal{O}}$ (note that this pseudo-Hamiltonian depends on both the operator and saddle-point state). Somewhat in parallel to the Quench Action itself, the so-defined operator pseudo-Hamiltonian is a measure of the importance of eigenstates as far as expectation values are concerned. Important states have small (real part of) pseudo-energy; conversely, states with negligible matrix elements correspond to high pseudo-energy states. The potentially crucial simplification comes if the weights can then be expressed (as done above for the differential overlaps) in terms of one-, two-,... body functions,

$$\begin{eqnarray}&&{{h}_{\mathbf{e}}}=\underset{i}{\sum}h\left({{e}_{i}}\right)+\underset{{{i}_{1}}<{{i}_{2}}}{\sum}{{h}^{(2)}}\left({{e}_{{{i}_{1}}}},{{e}_{{{i}_{2}}}}\right)+...\end{eqnarray} \tag{ 73 }$$

An operator which is not too dissimilar to the operator creating individual particle-hole excitations would naturally be represented using only the first few terms to achieve good accuracy. The advantage of looking for such a representation is that the whole operator structure is then encoded in the smallest number of parameters possible. This topic, which has not yet been implemented beyond unpublished isolated cases and goes beyond the Quench Action towards the computation of correlations on generic states, will be investigated further in future works.

To summarize, in order to reconstruct the full time dependence of the expectation value of the weak, smooth operator $\mathcal{O}$, the only ingredients that are needed are:

• (i)  
  the saddle-point distribution ${{\rho}_{\text{sp}}}$,
• (ii)  
  the excitation energy function⁶ ε
• (iii)  
  the characteristic quench overlap function s (and if needed its higher-body parts s⁽²⁾ etc), and
• (iv)  
  the matrix elements $\langle {{\rho}_{\text{sp}}}|\mathcal{O}|{{\rho}_{\text{sp}}};\left\{\mathbf{e}\right\}\rangle$ for states around the saddle point, either directly or in terms of the relevant operator's pseudo-Hamiltonian.

These can be viewed as the distilled 'effective parameters' encoding the whole time evolution of a weak, smooth operator after the quench. Note that the first three items are operator-independent.

It is worthwhile here to discuss how one should properly take the thermodynamic limit of all our constructions up to now. The point to remember is that the underlying integrable model provides a fully regularized theory at all energies, to all orders in system size, for all states and operators. The finite N equations thus hide all the necessary details of a regularization scheme which needs to be appended to the field theory in order to make it well-defined. Though we discuss this aspect here in the context of the Quench Action and its application, the same reasonings are applicable in much more generality.

As explained above, the actual calculation of (66) should be performed using a defined microscopic realization of the steady state, for example the 'maximally flat' one,

$$\begin{eqnarray}&&\underset{\mathbf{e}}{\sum}{{\text{e}}^{-\delta {{S}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]-\text{i}{{\omega}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]t}}\langle {{\rho}_{\text{sm}}}|\mathcal{O}|{{\rho}_{\text{sm}}};\mathbf{e}\rangle =\underset{\mathbf{e}}{\sum}{{\text{e}}^{-\delta {{S}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]-\text{i}{{\omega}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]t}}\langle \left\{{{\rho}_{i}}\right\};\left\{0\right\}|\mathcal{O}|\left\{{{\rho}_{i}}\right\};\left\{0\right\};\mathbf{e}\rangle.\end{eqnarray} \tag{ 74 }$$

Let us specify the structure of the sum over excitations a little bit more precisely. In our box regularization, let us as before use latin indices i to label the boxes. The difference between the bra and ket states in (74) consists in in-box 'entropy-like' modifications, accompanied by out-box 'dispersing' excitations taking the form of particles or holes drilled on the overall density profile $\left\{{{\rho}_{i}}\right\}$. It is convenient to separate these two types of excitations explicitly, and to effectively perform a summation over the former. Denoting the number of such dispersing particle (resp. hole) excitations as n_p (resp. n_h) and the box which they disperse to by ${{i}_{a}},a=1,...,{{n}_{p}}$ (resp. ${{\bar{i}}_{b}},b=1,...,{{n}_{h}}$), we can represent the sum over excitations in the example above as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\mathbf{e}}{\sum}{{\text{e}}^{-\delta {{S}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]-\text{i}{{\omega}_{\mathbf{e}}}\left[{{\rho}_{\text{sp}}}\right]t}}\langle \left\{{{\rho}_{i}}\right\};\left\{0\right\}|\mathcal{O}|\left\{{{\rho}_{i}}\right\};\left\{0\right\};\mathbf{e}\rangle \\ \,\,\,\,\,\,=\underset{{{n}_{p}}=0}{\overset{\infty}{\sum}}\underset{{{n}_{h}}=0}{\overset{\infty}{\sum}}\frac{1}{{{n}_{p}}!{{n}_{h}}!}\underset{{{\left\{{{i}_{a}}\right\}}_{{{n}_{p}}}};{{\left\{{{{\bar{i}}}_{b}}\right\}}_{{{n}_{h}}}}}{\sum}\left[\underset{a=1}{\overset{{{n}_{p}}}{\prod}}\underset{b=1}{\overset{{{n}_{h}}}{\prod}}(1-{{\delta}_{{{i}_{a}}{{{\bar{i}}}_{b}}}})\right]{{\text{e}}^{-{{{\sum}^{}}_{a}}\,{{s}_{{{i}_{a}}}}+{{{\sum}^{}}_{b}}\,{{s}_{{{{\bar{i}}}_{b}}}}-\text{i}\left({{{\sum}^{}}_{a}}\,{{\varepsilon}_{{{i}_{a}}}}-{{{\sum}^{}}_{b}}\,{{\varepsilon}_{{{{\bar{i}}}_{b}}}}\right)t}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\times \underset{\left\{{{c}_{i}}\right\}}{\sum}\langle \left\{{{\rho}_{i}}\right\};\left\{0\right\}|\mathcal{O}|\left\{{{\rho}_{i}}\right\}{\cup}^{}\left\{{{i}_{a}}\right\},\left\{{{{\bar{i}}}_{b}}\right\};\left\{{{c}_{i}}\right\}\rangle. \end{array}\end{eqnarray} \tag{ 75 }$$

A few things are worth clarifying. First of all, $\sum\nolimits_{\left\{{{c}_{i}}\right\}}$ represents as before the in-box summations at fixed density profile $\left\{{{\rho}_{i}}\right\}$. The factors $1-{{\delta}_{{{i}_{a}}{{{\bar{i}}}_{b}}}}$ enforce the convention that dispersing particle and hole excitations are chosen by convention not to land in the same box, this being meaningless since the density set $\left\{{{\rho}_{i}}\right\}$ is then not modified (the sum over configurations {c_i} already takes care of these terms). Multiple particles (resp. holes) can however land in the same box. We have also exploited the additivity of the differential overlap and energy functions around the saddle point.

It is convenient at this stage to introduce 'box-averaged' states

$$\begin{eqnarray}&&|\left\{{{\rho}_{i}}\right\}{{\rangle}_{b}}\equiv {{\text{e}}^{-\frac{1}{2}S_{\left\{{{\rho}_{i}}\right\}}^{\text{YY}}}}\underset{\left\{{{c}_{i}}\right\}}{\sum}|\left\{{{\rho}_{i}}\right\};\left\{{{c}_{i}}\right\}\rangle \end{eqnarray} \tag{ 76 }$$

obtained by uniformly summing over all in-box configurations at fixed density. Box-averaged states are not exact eigenstates of our Hamiltonian; nevertherless, their energy is well-defined with fluctuations which are exponentially small in system size and can be neglected in the thermodynamic limit. Box-averaged states are still orthonormal in the sense that

$$\begin{eqnarray}&&{{}_{b}}{{\langle \left\{\rho _{i}^{l}\right\}|\left\{\rho _{i}^{r}\right\}\rangle}_{b}}={{\delta}_{\left\{\rho _{i}^{l}\right\},\left\{\rho _{i}^{r}\right\}}}.\end{eqnarray} \tag{ 77 }$$

One can also introduce for convenience the usual ladder operators ${{z}_{j}},z_{j}^{\dagger}$ such that

$$\begin{eqnarray}&&z_{j}^{\dagger}|\left\{{{\rho}_{i}}\right\}{{\rangle}_{b}}=\,\left|\left\{{{\rho}_{i}}+\frac{{{\delta}_{ij}}}{{{l}_{j}}}\right\}\right\rangle_{b},\quad {{z}_{j}}|\left\{{{\rho}_{i}}\right\}{{\rangle}_{b}}=\,\left|\left\{{{\rho}_{i}}-\frac{{{\delta}_{ij}}}{{{l}_{j}}}\right\}\right\rangle_{b}\end{eqnarray} \tag{ 78 }$$

with simple algebra and additional lowest/highest-weight conditions

$$\begin{eqnarray}&&[{{z}_{j}},z_{{{j}^{\prime}}}^{\dagger}]={{\delta}_{j{{j}^{\prime}}}},\quad z_{j}^{\dagger}|\{{{\rho}_{i}}\}{{\rangle}_{b}}{{|}_{{{\rho}_{j}}={{l}_{j}}}}=0,\quad {{z}_{j}}|\{{{\rho}_{i}}\}{{\rangle}_{b}}{{|}_{{{\rho}_{j}}=0}}=0.\end{eqnarray} \tag{ 79 }$$

These simple operators are somewhat reminiscent of Zamolodchikov-Faddeev operators in integrable field theory [75] (these operators have indeed been used in the context of quenches, see e.g. [76, 77]). We however want to emphasize that they are not the same: the simple operators here are still fully microscopics-aware, and permit calculations in principle to all orders in system size (to put it in field theory language: there is no need for further regularization, there are no infinities present at this stage).

Introducing operators f_I, $f_{I}^{\dagger}$ defined as operators removing/adding an occupation at quantum number I in the full microscopic model, and the in-box symmetric combination ${{f}_{j}}\equiv \sum\nolimits_{I\in \text{box}~j}{{f}_{I}}$ (one can think of the (un)occupied quantum numbers as a pseudo-spin-1/2 degree of freedom, with operators ${{f}_{I}}/f_{I}^{\dagger}$ as a spin lowering/raising operator, and the ${{f}_{j}},f_{j}^{\dagger}$ as total spin operators in the representation with total spin l_j/2), we have the actions (taking the factors coming from the entropy into account)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} f_{j}^{\dagger}|\left\{{{\rho}_{i}}\right\}{{\rangle}_{b}}={{\left[\left({{n}_{j}}+1\right)\left({{l}_{n}}-{{n}_{j}}\right)\right]}^{1/2}}\left|\left\{{{\rho}_{i}}+\frac{{{\delta}_{ij}}}{{{l}_{j}}}\right\}\right\rangle_{b}, \\ {{f}_{j}}|\left\{{{\rho}_{i}}\right\}{{\rangle}_{b}}={{\left[{{n}_{j}}\left({{l}_{j}}-{{n}_{j}}+1\right)\right]}^{1/2}}\left|\left\{{{\rho}_{i}}-\frac{{{\delta}_{ij}}}{{{l}_{j}}}\right\}\right\rangle_{b} \end{array}\end{eqnarray} \tag{ 80 }$$

and thus the identities

$$\begin{eqnarray}&&z_{j}^{\dagger}=f_{j}^{\dagger}\frac{{{l}_{j}}}{{{\left[({{{\hat{\rho}}}_{j}}+\frac{1}{{{l}_{j}}})(1-{{{\hat{\rho}}}_{j}})\right]}^{1/2}}}~({{\rho}_{j}}\ne 1),\quad {{z}_{j}}={{f}_{j}}\frac{{{l}_{j}}}{{{\left[{{{\hat{\rho}}}_{j}}(1-{{{\hat{\rho}}}_{j}}+\frac{1}{{{l}_{j}}})\right]}^{1/2}}}~({{\rho}_{j}}\ne 0)\end{eqnarray} \tag{ 81 }$$

in which ${{\hat{\rho}}_{j}}$ returns the density in box j. The use of introducing box-averaged states is that these will become the well-defined states in the thermodynamic limit, parametrized by the density distribution. We can then define a box-renormalized operator ${{\mathcal{O}}^{\text{br}}}$, with matrix elements in the box-averaged state basis absorbing the renormalization coming from the in-box summation,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{1}{\prod\nolimits_{a=1}^{{{n}_{p}}}({{l}_{{{i}_{a}}}}-{{n}_{{{i}_{a}}}})\prod\nolimits_{b=1}^{{{n}_{h}}}{{n}_{{{{\bar{i}}}_{b}}}}}\underset{\left\{{{c}_{i}}\right\}}{\sum}\langle \left\{{{\rho}_{i}}\right\};\left\{0\right\}|\mathcal{O}|\left\{{{\rho}_{i}}\right\}{\cup}^{}\left\{{{i}_{a}}\right\},\left\{{{{\bar{i}}}_{b}}\right\};\left\{{{c}_{i}}\right\}\rangle \\ \,\,\,\,\,=\frac{{{\text{e}}^{-S_{\left\{{{\rho}_{i}}\right\}}^{\text{YY}}}}}{\prod\nolimits_{a=1}^{{{n}_{p}}}({{l}_{{{i}_{a}}}}-{{n}_{{{i}_{a}}}})\prod\nolimits_{b=1}^{{{n}_{h}}}{{n}_{{{{\bar{i}}}_{b}}}}}\underset{\left\{{{c}_{i}}\right\},\left\{c_{i}^{\prime}\right\}}{\sum}\langle \left\{{{\rho}_{i}}\right\};\left\{c_{i}^{\prime}\right\}|\mathcal{O}|\left\{{{\rho}_{i}}\right\}{\cup}^{}\left\{{{i}_{a}}\right\},\left\{{{{\bar{i}}}_{b}}\right\};\left\{{{c}_{i}}\right\}\rangle \\ \,\,\,\,\,\equiv {{}_{b}}{{\langle \left\{{{\rho}_{i}}\right\}|{{\mathcal{O}}^{\text{br}}}|\left\{{{\rho}_{i}}\right\};\left\{{{i}_{a}}\right\},\left\{{{{\bar{i}}}_{b}}\right\}\rangle}_{b}} \end{array}\end{eqnarray} \tag{ 82 }$$

in which we have in the first step exploited the assumption that operator $\mathcal{O}$ is smooth. The factors ${{l}_{i}}-{{n}_{i}}$ and ${{n}_{{\bar{i}}}}$ are included for convenience to compensate for the fact that a particle excitation can be positioned at ${{l}_{i}}-{{n}_{i}}$ positions in box i, and the hole in ${{n}_{{\bar{i}}}}$ positions, the microscopically-defined matrix element being assumed to be insensitive to the precise positioning of the dispersing excitations at the level of resolution corresponding to the box size (this holds true in general and can be verified on a case-by-base basis).

The point is that the renormalization factor is density profile-specific, but is expected to be only weakly operator dependent and does not depend on the number or positioning of the dispersing excitations (provided there is a denumerable number of those). The summation over in-box configurations is reminiscent of the summation over soft modes used in [78] (where it was done for soft modes around ground states, leading to an anomalous power-law renormalization of matrix elements; here, being performed around finite-entropy states, the renormalization can be of exponential size). Actually performing the summations representing the operator renormalization is a nontrivial task, which needs to be performed on a case-by-case basis. Note that the value of the matrix elements in this box-averaged state basis (right-hand side of (82)) is now not bounded in the thermodynamic limit. In particular, they can develop singularities when particle and hole excitations come closer together (namely when they disperse to boxes which are closer and closer to each other). Additionally, they also display a crossing symmetry, namely an equality between putting e.g. a hole in the bra or a corresponding particle in the ket. These matrix elements thus essentially obtain properties reminiscent of some of those of form factors in integrable field theory [79]. We will get back to this point in more detail in future publications.

Going back to our discussion, our summation (74) has now become

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{{{n}_{p}}=0}{\overset{\infty}{\sum}}\underset{{{n}_{h}}=0}{\overset{\infty}{\sum}}\frac{1}{{{n}_{p}}!{{n}_{h}}!}\underset{{{\left\{{{i}_{a}}\right\}}_{{{n}_{p}}}};{{\left\{{{{\bar{i}}}_{b}}\right\}}_{{{n}_{h}}}}}{\sum}\left[\underset{a=1}{\overset{{{n}_{p}}}{\prod}}\underset{b=1}{\overset{{{n}_{h}}}{\prod}}(1-{{\delta}_{{{i}_{a}}{{{\bar{i}}}_{b}}}})\right]\underset{a=1}{\overset{{{n}_{p}}}{\prod}}\left({{l}_{{{i}_{a}}}}-{{n}_{{{i}_{a}}}}\right)\underset{b=1}{\overset{{{n}_{h}}}{\prod}}{{n}_{{{{\bar{i}}}_{b}}}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times \,{{\text{e}}^{-{{{\sum}^{}}_{a}}{{s}_{{{i}_{a}}}}+{{{\sum}^{}}_{b}}{{s}_{{{{\bar{i}}}_{b}}}}-i\left({{{\sum}^{}}_{a}}{{\varepsilon}_{{{i}_{a}}}}-{{{\sum}^{}}_{b}}{{\varepsilon}_{{{{\bar{i}}}_{b}}}}\right)t}}{{}_{b}}{{\langle \left\{{{\rho}_{i}}\right\}|{{\mathcal{O}}^{\text{br}}}|\left\{{{\rho}_{i}}\right\};\left\{{{i}_{a}}\right\},\left\{{{{\bar{i}}}_{b}}\right\}\rangle}_{b}} \end{array}\end{eqnarray} \tag{ 83 }$$

in which

$$\begin{eqnarray}&&|\left\{{{\rho}_{i}}\right\};\left\{{{i}_{a}}\right\},\left\{{{\bar{i}}_{b}}\right\}{{\rangle}_{b}}=\underset{a}{\prod}z_{{{i}_{a}}}^{\dagger}\underset{b}{\prod}{{z}_{{{{\bar{i}}}_{b}}}}|\left\{{{\rho}_{i}}\right\}{{\rangle}_{b}}.\end{eqnarray} \tag{ 84 }$$

One can then proceed with the continuum limit, introduce integrals over the quantum numbers x  =  I/L by using $\Delta {{x}_{i}}=\frac{{{l}_{i}}}{L}$, namely $\sum\nolimits_{{{i}_{a}}}\left(...\right)={\int}^{}\text{d}x\frac{L}{{{l}_{i}}}\left(...\right)$, to get (in a self-explanatory notation)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{{{n}_{p}}=0}{\overset{\infty}{\sum}}\underset{{{n}_{h}}=0}{\overset{\infty}{\sum}}\frac{{{L}^{{{n}_{p}}+{{n}_{h}}}}}{{{n}_{p}}!{{n}_{h}}!}\fint \underset{a=1}{\overset{{{n}_{p}}}{\prod}}\text{d}{{x}_{a}}\underset{b=1}{\overset{{{n}_{h}}}{\prod}}\text{d}{{{\bar{x}}}_{b}}\underset{a=1}{\overset{{{n}_{p}}}{\prod}}{{\rho}^{h}}\left({{x}_{a}}\right)\underset{b=1}{\overset{{{n}_{h}}}{\prod}}\rho \left({{{\bar{x}}}_{b}}\right){{\text{e}}^{-\delta S\left(\mathbf{x},\mathbf{\bar{x}}\right)-\text{i}\omega \left(\mathbf{x},\mathbf{\bar{x}}\right)t}}\langle {{\rho}_{\text{sm}}}|{{\mathcal{O}}^{\text{br}}}|{{\rho}_{\text{sm}}};\mathbf{x},\mathbf{\bar{x}}\rangle \end{array}\end{eqnarray} \tag{ 85 }$$

in which $\fint$ means the integral with 'principal part' extraction of all points ${{x}_{a}}={{\bar{x}}_{b}}~\forall ~a,b$. Note that instead of integrating over quantum number space variables, one can equivalently integrate over rapidity variables by using the transformation function $x\left(\lambda \right)$ (17) and its Jacobian (18). As described earlier, the energy ω and differential overlap function $\delta S$ are also expressible in terms of sums over excitations, the energy being purely one-body, and the differential overlap perhaps having some higher-body parts. Note that formula (85) parallels equation (2.38) in [80], the difference resting in the precise definition of the matrix element (in our case here, box-averaged), and is an adaptation of the LeClair–Mussardo formula [81] to this highly-off-vacuum context.

Expression (85) and its mirror term as per (66) then give a more traditionally field theory-looking version of the fundamental equation of the Quench Action formalism, in the sense that the excitations one sums over are strictly denumerable. Though somewhat formal, they make clear which kind of contributions need to be taken into consideration when computing the time dependence in the thermodynamic limit. They might serve as a different way to regularize quenches in integrable field theory [50, 53, 63, 76, 77, 82–87].

The first problem treated [1] using the general formalism of the Quench Action was that of the transverse-field Ising model

$$\begin{eqnarray}&&H(h)=-J\underset{j=1}{\overset{L}{\sum}}[\sigma _{j}^{x}\sigma _{j+1}^{x}+h\sigma _{j}^{z}],\end{eqnarray} \tag{ 86 }$$

with J, h  >  0. This model has two phases, one with ferromagnetic order along the x direction for h  <  1, the other being paramagnetic with h  >  1, these being separated by a quantum critical point in the Ising universality class (see e.g. [88]). Bethe Ansatz is not needed here; diagonalization can be performed using Jordan-Wigner and Bogoliubov transformations, leading to free dispersive modes ${{\alpha}_{k}}$ with diagonal Hamiltonian

$$\begin{eqnarray}&&H(h)=\underset{k}{\sum}{{\varepsilon}_{n}}(k)\left(\alpha _{k}^{\dagger}{{\alpha}_{k}}-\frac{1}{2}\right),\quad {{\varepsilon}_{h}}(k)=2J\sqrt{1+{{h}^{2}}-2h\cos k}.\end{eqnarray} \tag{ 87 }$$

The first quench protocol considered consists in preparing the system in the ground state of (86) for an initial value h₀ in the paramagnetic phase h₀  >  1, and to suddenly switch the field to value h from t  =  0 onwards.

The post-quench reduced density matrix has been computed in [89], giving ${{\rho}^{\text{DM}}}(t)=|{{\rho}_{\text{sp}}}\rangle \langle {{\rho}_{\text{sp}}}|\propto {{\text{e}}^{\frac{1}{4}{{\mathbf{a}}^{T}}\mathbf{Wa}}}$ with $\mathbf{a}$ being a vector of Majorana modes ${{a}_{2n-1}}=\left[\prod\nolimits_{m<n}\sigma _{m}^{z}\right]\sigma _{n}^{x}$ , ${{a}_{2n}}=\left[\prod\nolimits_{m<n}\sigma _{m}^{z}\right]\sigma _{n}^{y}$ with anticommutation $\left\{{{a}_{m}},{{a}_{n}}\right\}=2{{\delta}_{mn}}$, the matrix W being given by $\tanh (W/2)=\Gamma$ [90, 91] in which

$$\begin{eqnarray}&&{{\Gamma}_{jk}}=~\text{Tr}~\left[{{\rho}^{\text{DM}}}(t){{a}_{k}}{{a}_{j}}\right]-{{\delta}_{jk}}.\end{eqnarray} \tag{ 88 }$$

The steady state can be obtained by applying the Quench Action method (we refer the reader to [1] for all details). In this case, the solution to the relevant GTBA can be obtained directly from the fact that the momentum occupations are conserved. This leads to the steady state distribution

$$\begin{eqnarray}&&\rho (k)=\frac{1-\cos {{\Delta}_{k}}}{4\pi},\quad \cos {{\Delta}_{k}}=4{{J}^{2}}\,\left(1+h{{h}_{0}}-\left(h+{{h}_{0}}\right)\cos k\right)/\left[{{\varepsilon}_{h}}(k){{\varepsilon}_{{{h}_{0}}}}(k)\right]\end{eqnarray} \tag{ 89 }$$

in which the momenta ${{\kappa}_{j}}$ are distributed according to the saddle-point density (89), ${{\kappa}_{j+1}}={{\kappa}_{j}}+\frac{1}{L\rho \left({{\kappa}_{j}}\right)}$. This leads to a Gaussian steady state $|{{\rho}_{\text{sp}}}\rangle =\prod\nolimits_{j}\alpha _{{{\kappa}_{j}}}^{\dagger}\alpha _{-{{\kappa}_{j}}}^{\dagger}|0;h\rangle$ (the state $|0;h\rangle$ being the vacuum at field value h) with density matrix $\rho _{\text{sp}}^{\text{DM}}=\,|{{\rho}_{\text{sp}}}\rangle \langle {{\rho}_{\text{sp}}}|$ completely characterized by a correlation matrix whose only nonzero elements are

$$\begin{eqnarray}&&{{({{\Gamma}_{\text{sp}}})}_{2l-1,2l-2n}}=-\frac{\text{i}}{L}\underset{k}{\sum}\frac{{{\text{e}}^{-\text{i}nk}}(h-{{\text{e}}^{\text{i}k}})(1-2{{\delta}_{k,{{\kappa}_{j}}}})}{\sqrt{1+{{h}^{2}}-2h\cos k}}.\end{eqnarray} \tag{ 90 }$$

This density matrix obeys ${{\Gamma}_{\text{sp}}}=\Gamma(t\to \infty )$ and thus indeed equation (67) holds for any (products of) local operators, offering a simple way of reproducing the form obtained beforehand [43–45]. For time-dependent expectation values, equation (66) can similarly be verified at the level of the (reduced) density matrix.

To summarize, although this quench problem had been explicitly solved before, the Quench Action logic allowed to reproduce steady-state results (and also to reobtain the known time-dependent relaxation behaviour) in a computationally less heavy fashion.

The previous example pertaining to the transverse-field Ising belongs to the class of problems which, due to the absence of interactions, are solvable without invoking the technology of the Bethe Ansatz. Although the dynamics of such models is surely representative of generic situations (in particular because the mapping between interesting observables and the diagonalized modes is sometimes very complex and nonlocal) and the GGE is expected to hold, it is interesting to study situations where interactions are inevitably strong and no mapping to a simple free theory exists.

Let us thus now turn our attention to what is to the best of the author's knowledge the first quench problem to be solved analytically in the thermodynamic limit (by which we mean obtaining an exact analytical characterization of the steady state), for post-quench time evolution in the presence of nontrivial interactions.

In the context of bosons in one dimension, if one takes as initial state the ground state of the noninteracting bosonic gas on a periodic interval of circumference L, namely the Bose–Einstein condensate-like⁷ state

$$\begin{eqnarray}&&|{{\Psi}_{0}}\rangle =\frac{1}{\sqrt{{{L}^{N}}N!}}{{(\psi _{k=0}^{\dagger})}^{N}}|0\rangle \end{eqnarray} \tag{ 91 }$$

in which $|0\rangle$ is the Fock vacuum and $\psi _{k=0}^{\dagger}$ creates a zero-momentum particle, one can obtain an interesting quench problem by simply turning interactions on from t  =  0 onwards. For definiteness, the interactions are taken to be ultralocal, i.e. the post-quench Hamiltonian is defined as that of the Lieb–Liniger gas [92],

$$\begin{eqnarray}&&{{H}_{\text{LL}}}=-\underset{j=1}{\overset{N}{\sum}}\frac{{{\partial}^{2}}}{\partial _{{{x}_{j}}}^{2}}+2c\underset{{{j}_{1}}<{{j}_{2}}}{\sum}\delta ({{x}_{{{j}_{1}}}}-{{x}_{{{j}_{2}}}})\end{eqnarray} \tag{ 92 }$$

with c  >  0. In this case, the Bethe equations take the form

$$\begin{eqnarray}&&{{\lambda}_{j}}+\frac{1}{L}\underset{l=1}{\overset{N}{\sum}}\phi ({{\lambda}_{j}}-{{\lambda}_{l}})=\frac{2\pi}{L}{{I}_{j}},\quad j=1,...,N,\end{eqnarray} \tag{ 93 }$$

where $\phi (\lambda )\equiv 2~\text{atan}(\lambda /c)$ and the quantum numbers are distinct half-odd integers for N even and integers for N odd. The wavefunction itself is given by the Bethe Ansatz,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\Psi}_{N}}(\{x\}|\{\lambda \}) & = & \underset{N\geqslant {{j}_{1}}>{{j}_{2}}\geqslant 1}{\prod}~\text{sgn}({{x}_{{{j}_{1}}}}-{{x}_{{{j}_{2}}}})~\text{sgn}({{\lambda}_{{{j}_{1}}}}-{{\lambda}_{{{j}_{2}}}}) \\ {} & {} & \times \underset{P\in {{\pi}_{N}}}{\sum}{{(-1)}^{[P]}}{{\text{e}}^{\text{i}\sum\nolimits_{j=1}^{N}{{\lambda}_{{{P}_{j}}}}{{x}_{j}}+\frac{\text{i}}{2}{{{\sum}^{}}_{N\geqslant {{j}_{1}}>{{j}_{2}}\geqslant 1}}~\text{sgn}({{x}_{{{j}_{1}}}}-{{x}_{{{j}_{2}}}})\phi ({{\lambda}_{{{P}_{{{j}_{1}}}}}}-{{\lambda}_{{{P}_{{{j}_{2}}}}}})}}. \end{array}\end{eqnarray} \tag{ 94 }$$

In the initial BEC state, the local number density fluctuations are large, and the momentum distribution function is by definition a delta function at zero momentum carrying the weight of all particles; turning repulsive interactions on must then lead among others to a broadening of the momentum distribution function, with the suddenly quenched interaction energy being partly converted to kinetic energy as time progresses. Since all energies are incommensurately related at finite c, one does not expect any persistent oscillations or recurrences in the thermodynamic limit, but rather to see the system effectively relax to a steady state due to quantum dephasing.

The search for a solution to this quench problem has an interesting history (which will be discussed later on). It was solved exactly in [36] using the Quench Action formalism, and what follows is a brief summary of some important aspects of this solution.

4.2.1. Overlaps

In order to implement the Quench Action for this protocol, the overlaps between the BEC state (91) and the eigenstates (94) of the Lieb–Liniger gas (92) should be obtained. This seems, a priori, to pose insurmountable difficulties: although Bethe states are indeed composed of linear combinations of plane waves, which seem to correspond seamlessly to the eigenstates of a free Hamiltonian (and even more to the especially simple and structureless BEC state), the problem is that these plane waves are defined only on finite segments in coordinate space, and the overlap integrals always leave incommensurate phases hanging. A basic attempt at calculating the overlap using real-space integrals thus yields a factorially large sum of terms with no immediately obvious simplification in the general case. For the Lieb–Liniger gas in the Tonks–Girardeau limit, overlaps with the BEC state can be obtained [93] as products of the inverse of the rapidities.

Amazingly, for the BEC state, at any value of the interaction parameter, a dramatic simplification in fact turns out to be possible: one finds [36, 94] that the exact overlaps (which are only nonvanishing for parity-symmetric states, namely states in which the rapidity distribution is mirror-symmetric about the origin; this fact, with the Bethe equations, allows to perform simplifications) are exactly given by the relatively simple expression

$$\begin{eqnarray}&&\langle {{\Psi}_{0}}|{{\left\{\lambda \right\}}_{N/2}}{\cup}^{}{{\left\{-\lambda \right\}}_{N/2}}\rangle ={{\left[\frac{{{(cL)}^{-N}}N!}{\text{de}{{\text{t}}_{N}}{{G}_{jk}}}\right]}^{1/2}}\frac{\text{de}{{\text{t}}_{N/2}}G_{jk}^{Q}}{\prod\nolimits_{j=1}^{N/2}\frac{{{\lambda}_{j}}}{c}{{\left[\frac{\lambda _{j}^{2}}{{{c}^{2}}}+\frac{1}{4}\right]}^{1/2}}}\end{eqnarray} \tag{ 95 }$$

in which G_jk is the Gaudin matrix

$$\begin{eqnarray}&&{{G}_{jk}}={{\delta}_{jk}}\left[L+\underset{l=1}{\overset{N}{\sum}}K\left({{\lambda}_{j}},{{\lambda}_{l}}\right)\right]-K\left({{\lambda}_{j}},{{\lambda}_{k}}\right)\end{eqnarray} \tag{ 96 }$$

with kernel $K\left(\lambda,{{\lambda}^{\prime}}\right)=\frac{2c}{{{\left(\lambda -{{\lambda}^{\prime}}\right)}^{2}}+{{c}^{2}}}$, and G^Q is the Gaudin-like matrix

$$\begin{eqnarray}&&G_{jk}^{Q}={{\delta}_{jk}}\left[L+\underset{l=1}{\overset{N/2}{\sum}}{{K}^{Q}}\left({{\lambda}_{j}},{{\lambda}_{l}}\right)\right]-{{K}^{Q}}\left({{\lambda}_{j}},{{\lambda}_{k}}\right)\end{eqnarray} \tag{ 97 }$$

with ${{K}^{Q}}\left(\lambda,{{\lambda}^{\prime}}\right)=K\left(\lambda,{{\lambda}^{\prime}}\right)+K\left(\lambda,-{{\lambda}^{\prime}}\right)$. This overlap formula is extremely economical, in the sense that its complexity scales only with the third power of system size due to its determinant structure, instead of the factorially large worst-case expectation. It is reminiscent not only of the Gaudin norm formula, but of the matrix elements of local density or field operators for this model [26, 59, 61, 67].

4.2.2. The steady state from the Quench Action

For the application of the Quench Action, only the extensive part of the (logarithmic) overlap is required, which can be extracted from (95). The distribution-dependent extensive part of the overlap is

$$\begin{eqnarray}&&{{S}^{Q}}/L={\int}_{0}^{\infty}\text{d}\lambda \rho \left(\lambda \right)\log \,\left(\frac{{{\lambda}^{2}}}{{{c}^{2}}}\left(\frac{1}{4}+\frac{{{\lambda}^{2}}}{{{c}^{2}}}\right)\right).\end{eqnarray} \tag{ 98 }$$

This is then combined with the Yang–Yang entropy

$$\begin{eqnarray}&&{{S}^{\text{YY}}}/L={\int}_{0}^{\infty}\text{d}\lambda \left[{{\rho}^{t}}\left(\lambda \right)\log {{\rho}^{t}}\left(\lambda \right)-\rho \left(\lambda \right)\log \rho \left(\lambda \right)-{{\rho}^{h}}\left(\lambda \right)\log {{\rho}^{h}}\left(\lambda \right)\right]\end{eqnarray} \tag{ 99 }$$

(in which ${{\rho}^{t}}=\rho +{{\rho}^{h}}$; note that this expression for the entropy reflects the required parity-invariance of the states by integrating over the positive rapidity half-line only), to give the Quench Action (50). The variational condition (52) then gives (with proper handling of the filling of the gas, see the original paper [36] for details) a functional equation for the steady-state root distribution ${{\rho}_{\text{sp}}}$. Defining $\eta \left(\lambda \right)=\frac{{{\rho}^{h}}\left(\lambda \right)}{\rho \left(\lambda \right)}$, this saddle-point equation reads

$$\begin{eqnarray}&&\ln \eta \left(\lambda \right)=g\left(\lambda \right)-h-{\int}_{-\infty}^{\infty}\frac{\text{d}{{\lambda}^{\prime}}}{2\pi}K\left(\lambda -{{\lambda}^{\prime}}\right)\ln \left[1+{{\eta}^{-1}}\left({{\lambda}^{\prime}}\right)\right]\end{eqnarray} \tag{ 100 }$$

and is thus of GTBA form with driving term

$$\begin{eqnarray}&&g\left(\lambda \right)=\ln \left[\frac{{{\lambda}^{2}}}{{{c}^{2}}}\left(\frac{{{\lambda}^{2}}}{{{c}^{2}}}+\frac{1}{4}\right)\right]\end{eqnarray} \tag{ 101 }$$

and chemical potential h adjusted to satisfy the filling condition ${\int}^{}\text{d}\lambda \rho \left(\lambda \right)=n$. This GTBA equation can amazingly be solved analytically in terms of modified Bessel functions of the first kind,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\rho}_{\text{sp}}}\left(\lambda \right)=-\frac{\gamma}{4\pi}\frac{1}{1+a\left(\lambda \right)}\frac{\partial a\left(\lambda \right)}{\partial \gamma},\quad a\left(\lambda \right)\equiv \frac{2\pi}{\frac{\lambda}{n}\sinh \frac{2\pi \lambda}{c}}\,{{I}_{1-2\text{i}\frac{\lambda}{c}}}\,\left(\frac{4}{\sqrt{\gamma}}\right)\,{{I}_{1+2\text{i}\frac{\lambda}{c}}}\,\left(\frac{4}{\sqrt{\gamma}}\right). \end{array}\end{eqnarray} \tag{ 102 }$$

As expected, the plot of these distributions (see figure 1 of [36]) shows an origin-centered, broadened peak which is non-thermal in shape. From this distribution, steady-state properties can be computed. In [36], the local static density moments g_2,3 with

$$\begin{eqnarray}&&{{g}_{l}}=\langle {{\rho}_{\text{sp}}}|{{\left[{{\psi}^{\dagger}}(0)\right]}^{l}}\psi {{(0)}^{l}}|{{\rho}_{\text{sp}}}\rangle /{{n}^{l}}\end{eqnarray} \tag{ 103 }$$

were computed in the thermodynamic limit by using the above-obtained steady state and the method described in [63, 95]. The static structure factor

$$\begin{eqnarray}&&S(x)=\langle {{\rho}_{\text{sp}}}|\rho (x)\rho (0)|{{\rho}_{\text{sp}}}\rangle \end{eqnarray} \tag{ 104 }$$

was itself computed in a finite-size regularization of the steady state by using ABACUS [96] and an adaptation of the finite-temperature method developed in [97].

A direct, coordinate Bethe Ansatz-based numerical verification of the Quench Action predictions for the properties of the steady state was obtained in [98, 99]. Few-particle properties were well reproduced, 3-body ones less so, which is to be expected in view of finite-size corrections (the numerics was done for up to 5 particles; a 3-body term is then not a weak operator).

4.2.3. Charges, divergences and GGE

Let us end this section with a reminder of the interesting history of the (attempts at a) solution of this quench protocol. The natural starting point, namely to attempt a GGE treatment, was pursued in [100] based on the second-quantized representations of the conserved charges, the idea being that these can in principle be readily evaluated in the initial BEC state. For the Lieb–Liniger model, the natural charges to consider are those coming from the transfer matrix with evaluation point around infinity (see [21] and references therein). This leads to a hierarchy of charges ${{\hat{Q}}^{(n)}}$ with eigenvalue equations

$$\begin{eqnarray}&&{{\hat{Q}}^{(n)}}|{{\left\{\lambda \right\}}_{N}}\rangle ={{Q}^{(n)}}\left({{\left\{\lambda \right\}}_{N}}\right),\quad {{Q}^{(n)}}\left({{\left\{\lambda \right\}}_{N}}\right)=\underset{j=1}{\overset{N}{\sum}}\lambda _{j}^{n}.\end{eqnarray} \tag{ 105 }$$

The generic form for the GGE would then be, in operator and eigenvalue forms,

$$\begin{eqnarray}&&\underset{n=0}{\overset{\infty}{\sum}}{{\beta}_{n}}{{\hat{Q}}^{(n)}},\quad L\underset{n}{\sum}{{\beta}_{n}}{\int}^{}\text{d}\lambda \rho \left(\lambda \right){{\lambda}^{n}},\end{eqnarray} \tag{ 106 }$$

with parameters ${{\beta}_{n}}$ determined from the initial conditions as per the GGE prescription.

The implementation of the GGE for this problem however hits a snag. Higher conserved charges [101, 102] have terms which are not quantum mechanically normal-ordered. Although evaluating these charges (defined with c -dependent coefficients) on an eigenstate of the 'correct' Hamiltonian (namely the same c as that used in the definition of the charge) leads to a well-defined eigenvalue, due to mutual cancellations between delta functions originating from the kinetic and interaction terms. This cancellation however does not occur when evaluating a charge defined at ${{c}^{\prime}}$ on a state defined at $c\ne {{c}^{\prime}}$.

This is reflected in the divergence of the t  =  0 expectation value of the higher charges on the BEC state, as can be seen from direct calculations [103]. The exact Quench Action solution for the steady state reproduces these divergences correctly: for all post-quench interaction values c  >  0, the saddle-point distribution (102) has a $1/{{\lambda}^{4}}$ tail at large rapidity,

$$\begin{eqnarray}&&{{\rho}_{\text{sp}}}\left(\lambda \right)=\frac{1}{2\pi}\left[\frac{{{n}^{4}}{{\gamma}^{2}}}{{{\lambda}^{4}}}+\frac{{{n}^{6}}{{\gamma}^{3}}\left(24-\gamma \right)}{4{{\lambda}^{6}}}+...\right]\end{eqnarray} \tag{ 107 }$$

(these first two terms in the expansion were obtained in [104] using q-bosons; the Quench Action solution can be expanded to any order if one feels motivated). The presence of this tail has dramatic consequences as far as conserved charges are concerned, since now all even charges with index $\geqslant 4$ become infinite, ${\int}^{}{{\rho}_{\text{sp}}}\left(\lambda \right){{\lambda}^{2n}}\to \infty$, $n\geqslant 2$. The divergence of Q⁽⁴⁾ is like $\delta (x=0)$ (namely: the dimension of momentum space, infinite here since we are in the continuum without a UV cutoff), with higher charges diverging more and more strongly.

The GGE logic therefore cannot be applied to this particular quench problem (which is another reason, at least for the author, to find this problem particularly interesting). These presence of these divergences motivated the study in [104] (updating [100]) by considering a regularization of the problem in terms of q-bosons, for which a partial GGE could be consistently implemented. A full GGE was however out of reach, and it remained for the Quench Action to provide for a full solution.

This problem with divergences is not confined to the BEC to Lieb–Liniger quench. The three preconditions for its existence are that the Hilbert space be infinite-dimensional (e.g. here on the continuum interval, there is no UV cutoff), that the maximal value of the interaction potential be unbounded (here, a delta function), and that the quench protocol involve a change of this interaction potential. It is thus anticipated that a more generic Lieb–Liniger quench protocol going from c to a different ${{c}^{\prime}}$ would display the same divergences, and also any form in interaction quenches in multi-species generalizations of Lieb–Liniger.

Using the perspective offered by the exact solution of the problem using the Quench Action, we can note that the forms (106) and (98) are incompatible due to the logarithmic singularity in the latter (this being an exact result), which cannot be recaptured by power-sum (or polynomial-type) conserved charge eigenvalues. One suggestive way to 'repair' the GGE for the BEC to Lieb–Liniger quench is thus to include a charge which does not directly originate from the usual trace identities coming from the transfer matrix. If one included a 'log' charge with eigenvalue

$$\begin{eqnarray}&&{{Q}^{\text{log}}}(\{{{\lambda}_{j}}\})=\underset{j=1}{\overset{N}{\sum}}\ln [\lambda _{j}^{2}(\lambda _{j}^{2}+{{(c/2)}^{2}})],\end{eqnarray} \tag{ 108 }$$

then the exact QA free energy (98) trivially becomes of GGE form (with only this single charge!). This charge (a kind of 'logarithm of particle energy') does not appear physically meaningful at first sight, though it is mathematically completely well defined (its matrix elements on the basis of Bethe states are all well-defined, at least for the even particle number sector) and it gives meaningful, extensive values on 'normal-looking' Bethe states and is thus probably of (quasi?)local nature. This way of 'repairing' the GGE does not correspond to a kind of UV regularization as offered by the q -boson approach. Whether one can ascribe any meaning to all of this is an open issue.

A physically perhaps more intuitive way to save the situation is to consider a slightly modified quench problem, whereby the initial state is not the perfect BEC state, but a UV/high-energy regularized one ${{\text{e}}^{-\epsilon {{H}_{\text{LL}}}}}|{{\Psi}_{0}}\rangle$. This modification of the initial state leads to the exponential suppression of the mode occupation at rapidities $\lambda \gg \frac{1}{\sqrt{\epsilon}}$, and thus to the regularization of all infinities in the initial values of the charges. This approach, which is also a kind of UV regularization, is discussed in more detail in [105].

4.2.4. Time evolution

The time evolution following the BEC to Lieb–Liniger quench can in principle be computed using the Quench Action method, using the fundamental representation (66). Besides the saddle-point distribution (102), one needs the differential overlap function (64), which was calculated in the original paper [36] (see equation (A14)) and indeed obeys the expected decomposition (69) as a sum of one-body terms,

$$\begin{eqnarray}&&\delta {{S}_{\mathbf{e}}}[{{\rho}_{\text{sp}}}]=\underset{k=1}{\overset{n}{\sum}}[\delta s(\tilde{{\lambda}}_{k}^{p})-\delta s(\tilde{{\lambda}}_{k}^{h})],\end{eqnarray} \tag{ 109 }$$

where ${{\tilde{\lambda}}^{p,h}}$ are the particle and hole excitation rapidities and

$$\begin{eqnarray}&&\delta s\left({{\tilde{{\lambda}}}^{p}}\right)-\delta s\left({{\tilde{{\lambda}}}^{h}}\right)={\int}_{0}^{\infty}\text{d}\lambda {{\rho}_{\text{sp}}}\left(\lambda \right)\frac{1+8\frac{{{\lambda}^{2}}}{{{c}^{2}}}}{\lambda \left(1+4\frac{{{\lambda}^{2}}}{{{c}^{2}}}\right)}F(\lambda |{{\tilde{{\lambda}}}^{p}},{{\tilde{{\lambda}}}^{h}})+\ln \left(\frac{{{\tilde{{\lambda}}}^{p}}\sqrt{{{\left({{\tilde{{\lambda}}}^{p}}/c\right)}^{2}}+1/4}}{{{\tilde{{\lambda}}}^{h}}\sqrt{{{\left({{\tilde{{\lambda}}}^{h}}/c\right)}^{2}}+1/4}}\right)\end{eqnarray} \tag{ 110 }$$

with backflow function obeying

$$\begin{eqnarray}&&2\pi F\left(\lambda |{{\tilde{\lambda}}^{p}},{{\tilde{\lambda}}^{h}}\right)-{\int}_{-\infty}^{\infty}\text{d}{{\lambda}^{\prime}}K\left(\lambda -{{\lambda}^{\prime}}\right){{\vartheta}_{\text{sp}}}\left({{\lambda}^{\prime}}\right)F\left({{\lambda}^{\prime}}|{{\tilde{\lambda}}^{p}},{{\tilde{\lambda}}^{h}}\right)=\phi \left(\lambda -{{\tilde{\lambda}}^{p}}\right)-\phi \left(\lambda -{{\tilde{\lambda}}^{h}}\right)\end{eqnarray} \tag{ 111 }$$

in which ${{\vartheta}_{\text{sp}}}\left(\lambda \right)={{\left[1+\rho _{\text{sp}}^{h}\left(\lambda \right)/{{\rho}_{\text{sp}}}\left(\lambda \right)\right]}^{-1}}$ is the saddle-point state's filling function.

These ingredients were used in the specific case of a quench from the BEC state to hard-core (Tonks–Girardeau) bosons, to obtain the exact analytical expression for the time-evolved one-body density matrix [106], this taking the form of a difference between two Fredholm Pfaffians. In addition, the previously-derived time-dependent density-density correlation [54] was reobtained using the Quench Action logic.

The QA was further employed to study the time dependence of observables at generic interaction in [107], suggesting that power-law relaxation should be generic at late times for observables which do not create large numbers of particle-hole excitations. Further investigation of specific examples of other observables is a promising research direction for the future.

Following up on the solution to the BEC to Lieb–Liniger quench described above, the Quench Action was then applied to a problem in the context of spin chains. More precisely, the initial state defined by the Néel state

$$\begin{eqnarray}&&|{{\Psi}_{0}}\rangle \equiv \frac{1}{\sqrt{2}}\left(|\uparrow \downarrow \uparrow \downarrow...\rangle \,+\,|\downarrow \uparrow \downarrow \uparrow...\rangle \right)\end{eqnarray} \tag{ 112 }$$

was time-evolved using the XXZ Hamiltonian

$$\begin{eqnarray}&&H=J\underset{j=1}{\overset{N}{\sum}}[S_{j}^{x}S_{j+1}^{x}+S_{j}^{y}S_{j+1}^{y}+\Delta S_{j}^{z}S_{j+1}^{z}]\end{eqnarray} \tag{ 113 }$$

with anisotropy parameter $\Delta \geqslant 1$ lying in the antiferromagnetic region. The initial state (112) is the ground state of the Ising Hamiltonian at $\Delta \to \infty$, so this can be viewed as an 'interaction' quench in the spin chain language. For definiteness, periodic boundary conditions are used.

In contrast to the repulsive Lieb–Liniger model, the XXZ chain admits bound states, taking the form of string patterns in the solution to the Bethe equations,

$$\begin{eqnarray}&&\lambda _{\alpha}^{j,a}=\lambda _{\alpha}^{j}+\text{i}\frac{\eta}{2}(n+1-2a)+\text{i}\delta _{\alpha}^{j,a},\quad a=1,...,n,\end{eqnarray} \tag{ 114 }$$

η being related to the chain's anisotropy parameter by the relation $\Delta =\cosh \eta$. The real parameter $\lambda _{\alpha}^{j}$ represents the string center, namely the 'center of mass' of the composite object represented by the n rapidities. The string deviations $\delta _{\alpha}^{j,a}$ are in most circumstances exponentially small in system size, and can typically be neglected, in which case one works in the so-called string hypothesis. Under this, the Bethe equations become equations for the string centers, the Bethe–Gaudin–Takahashi (BGT) equations [20, 22, 108, 109]

$$\begin{eqnarray}&&{{\theta}_{j}}(\lambda _{\alpha}^{j})-\frac{1}{N}\underset{k=1}{\overset{{{N}_{s}}}{\sum}}\underset{\beta =1}{\overset{{{M}_{k}}}{\sum}}{{\Theta}_{jk}}(\lambda _{\alpha}^{j}-\lambda _{\beta}^{k})=\frac{2\pi}{N}I_{\alpha}^{j}\end{eqnarray} \tag{ 115 }$$

in which

$$\begin{eqnarray}&&{{\theta}_{j}}\left(\lambda \right)=2\text{atan}\,\frac{\tan \left(\lambda \right)}{\tanh \left(n\eta /2\right)}+2\pi \left\lfloor \frac{\lambda}{\pi}+\frac{1}{2}\right\rfloor \end{eqnarray} \tag{ 116 }$$

and the string–string scattering phase shift is

$$\begin{eqnarray}&&{{\Theta}_{jk}}\left(\lambda \right)=\left(1-{{\delta}_{jk}}\right){{\theta}_{|\,j-k|}}\,\left(\lambda \right)+2{{\theta}_{|\,j-k|+2}}\,\left(\lambda \right)+...+2{{\theta}_{j+k-2}}\,\left(\lambda \right)+{{\theta}_{j+k}}\,\left(\lambda \right).\end{eqnarray} \tag{ 117 }$$

Unlike in the Lieb–Liniger model (where there is no UV cutoff and thus no limit on the quantum numbers), there exist limiting quantum numbers at each string level, which depend on the fillings of each level (we omit those details here for the sake of brevity). The string center rapidities are defined in the interval $[-\frac{\pi}{2},\frac{\pi}{2}]$.

The thermodynamic limit of the XXZ chain is expressed as before in terms of densities, this time each string type having its own density function. For example, the expressions for the energy and magnetization are

$$\begin{eqnarray}&&\frac{E}{N}=-\frac{h}{2}+\underset{n=1}{\overset{\infty}{\sum}}{\int}_{-\infty}^{\infty}\text{d}\lambda ~{{d}_{n}}(\lambda ){{\rho}_{n}}(\lambda ),\quad \frac{S_{tot}^{z}}{N}=\frac{1}{2}-\underset{n=1}{\overset{\infty}{\sum}}{\int}_{-\infty}^{\infty}\text{d}\lambda ~n{{\rho}_{n}}(\lambda )\end{eqnarray} \tag{ 118 }$$

where the thermal driving term is ${{d}_{n}}(\lambda )=\beta (hn-\pi J{{a}_{n}}(\lambda ))$, ${{a}_{n}}(\lambda )=\frac{1}{2\pi}\frac{\text{d}}{\text{d}\lambda}{{\theta}_{n}}(\lambda )$, while the expression for the Yang–Yang entropy is

$$\begin{eqnarray}&&{{S}^{\text{YY}}}[\{{{\rho}_{n}}\}]/N=\underset{n=1}{\overset{\infty}{\sum}}{\int}_{-\pi /2}^{\pi /2}\text{d}\lambda [\rho _{n}^{t}\ln \rho _{n}^{t}-{{\rho}_{n}}\ln {{\rho}_{n}}-\rho _{n}^{h}\ln \rho _{n}^{h}].\end{eqnarray} \tag{ 119 }$$

The usual TBA thermodynamic equilibrium equations, obtained by the saddle-point evaluation of the partition function's functional integral representation, take the form

$$\begin{eqnarray}&&\ln {{\eta}_{n}}(\lambda )=\beta {{d}_{n}}(\lambda )+\underset{m=1}{\overset{\infty}{\sum}}{{A}_{nm}}\ast \ln (1+\eta _{m}^{-1}(\lambda ))\end{eqnarray} \tag{ 120 }$$

in which ${{A}_{nm}}(\lambda )=\frac{1}{2\pi}\frac{\text{d}}{\text{d}\lambda}{{\Theta}_{nm}}(\lambda )$ and ${{\eta}_{n}}(\lambda )\equiv \rho _{n}^{h}(\lambda )/{{\rho}_{n}}(\lambda )$.

4.3.1. Overlaps

As before, the necessary starting point when aiming to implement a Quench Action treatment is to study the overlaps between the initial state (112) and the exact eigenstates of (113). Building on earlier work [110, 111] (see also [112]), it was shown in [113, 114] that these overlaps are given by [113]

$$\begin{eqnarray}&&\frac{\langle {{\Psi}_{0}}|\left\{\pm {{\lambda}_{j}}\right\}_{j=1}^{M/2}\rangle}{\Vert \left\{\pm {{\lambda}_{j}}\right\}_{j=1}^{M/2}\Vert}=\sqrt{2}\left[\underset{j=1}{\overset{M/2}{\prod}}\frac{\sqrt{\tan \left({{\lambda}_{j}}+\text{i}\eta /2\right)\tan \left({{\lambda}_{j}}-\text{i}\eta /2\right)}}{2\sin \left(2{{\lambda}_{j}}\right)}\right]\sqrt{\frac{\det\nolimits_{M/2}\left(G_{jk}^{+}\right)}{\det\nolimits_{M/2}\left(G_{jk}^{-}\right)}}\end{eqnarray} \tag{ 121 }$$

where

$$\begin{eqnarray}&&G_{jk}^{\pm}={{\delta}_{jk}}\left(N{{K}_{\eta /2}}\left({{\lambda}_{j}}\right)-\underset{l=1}{\overset{M/2}{\sum}}K_{\eta}^{+}\left({{\lambda}_{j}},{{\lambda}_{l}}\right)\right)+K_{\eta}^{\pm}\left({{\lambda}_{j}},{{\lambda}_{k}}\right),\quad j,k=1,\ldots,M/2\end{eqnarray} \tag{ 122 }$$

in which $K_{\eta}^{\pm}\left(\lambda,\mu \right)={{K}_{\eta}}\left(\lambda -\mu \right)\pm {{K}_{\eta}}\left(\lambda +\mu \right)$, and ${{K}_{\eta}}\left(\lambda \right)=\frac{\sinh \left(2\eta \right)}{\sin \left(\lambda +\text{i}\eta \right)\sin \left(\lambda -\text{i}\eta \right)}$. Note that a 'doubled' Gaudin matrix once again appears, this being reminiscent of the BEC overlaps (95). This is natural in view of the fact that the latter can in fact be viewed as a limit of the former (see [94]). Also, once again, only parity-invariant states contribute.

4.3.2. The steady state from the Quench Action

These ingredients were used in [115] to implement the Quench Action protocol and obtain an exact solution to the Néel to XXZ quench problem (an extended treatment is published in [116]). Let us recall a few of the important aspects of this solution.

The extensive part of the overlaps are extracted from (121), verifying that string deviations do not ruin the limit. These give us the overlap functional (46) as

$$\begin{eqnarray}&&{{S}^{o}}\left[\left\{\rho \right\}\right]=-\frac{N}{2}\underset{n=1}{\overset{\infty}{\sum}}{\int}_{0}^{\pi /2}\text{d}\lambda {{\rho}_{n}}\left(\lambda \right)\ln {{W}_{n}}\left(\lambda \right)\end{eqnarray} \tag{ 123 }$$

with ${{W}_{n}}\left(\lambda \right)$ being rapidity- and anisotropy-dependent functions explicitly given in [115], equations (S7) and (S8).

Everything is then in position to enforce the saddle-point condition (52), now generalized to many particle (string) types labeled by index n,

$$\begin{eqnarray}&&\left. \frac{\delta {{S}^{Q}}\left[\left\{\rho \right\}\right]}{\partial {{\rho}_{n}}\left(\lambda \right)}\right|_{\left\{{{\rho}_{n}}\right\}=\left\{{{\rho}_{n,\text{sp}}}\right\}}=0,\end{eqnarray} \tag{ 124 }$$

yielding [115] equations similar to the thermal equilibrium ones (120) but with different driving terms (h being again a chemical potential/magnetic field enforcing the filling constraint of zero magnetization):

$$\begin{eqnarray}&&\ln {{\eta}_{n}}\left(\lambda \right)=2n\left(\ln 4-h\right)+{{g}_{n}}\left(\lambda \right)+\underset{m=1}{\overset{\infty}{\sum}}{{A}_{nm}}\ast \ln \left(1+\eta _{m}^{-1}\left(\lambda \right)\right),\end{eqnarray} \tag{ 125 }$$

with drivings again given by logarithms,

$$\begin{eqnarray}&&{{g}_{n}}\left(\lambda \right)=\underset{l=0}{\overset{n-1}{\sum}}\ln \left[\frac{{{\sin}^{2}}\,\left(2\lambda \right)+{{\sinh}^{2}}\,\left[\eta \left(n-1-2l\right)\right]}{4|\tan \,\left[\lambda +\text{i}\frac{\eta}{2}(n-2l)\right]{{|}^{2}}}\right].\end{eqnarray} \tag{ 126 }$$

Similarly to the thermal equilibrium case, there exists also a more aesthetic form of these coupled equations making use of Takahashi's decoupling scheme. We omit these here for succinctness and refer to [115, 116] for this alternative version.

4.3.3. GGE using the transfer matrix charges

A GGE treatment including the traditional charges (those emanating from the usual 'spin-1/2' transfer matrix of the form (22), quasilocal charges not being known at the time) was proposed independently in [117] and [118] (the results of the latter subsequently extended in [119]). The GGE was thereafter exactly implemented using all these charges in [115]. A rather surprising fact was obtained in [115, 116], namely that the constraints

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}\langle \Psi(t)|{{\hat{Q}}^{(a)}}|\Psi(t)\rangle =\langle \left\{{{\rho}_{n,\text{sp}}}\right\}|{{\hat{Q}}^{(a)}}|\left\{{{\rho}_{n,\text{sp}}}\right\}\rangle \end{eqnarray} \tag{ 127 }$$

do not fix all the densities of the steady state, but rather only fix the density of holes of one-strings,

$$\begin{eqnarray}&&\rho _{1;sp}^{h}\left(\lambda \right)=\frac{{{\pi}^{2}}a_{1}^{3}\,\left(\lambda \right){{\sin}^{2}}\,\left(2\lambda \right)}{{{\pi}^{2}}a_{1}^{2}\,\left(\lambda \right){{\sin}^{2}}\left(2\lambda \right)+{{\cosh}^{2}}\,\left(\eta \right)}.\end{eqnarray} \tag{ 128 }$$

This means that all other densities corresponding to higher strings are actually left free to fluctuate, meaning that constraints seem to be missing, as also noticed in [120, 121].

This fact is indeed demonstrated by comparison to the exact Quench Action solution to the problem, leading to the somewhat shocking and unexpected conclusion that the GGE based on these traditional charges only, failed to converge to the exact Quench Action prediction for the steady state, the difference being markedly illustrated in the different distributions of even-length bound states at the origin (see figure 1 in [115]). These differences in distributions affect the observables, for example the local spin correlations; these differences were verified using numerics based on a linked-cluster expansion. In a back-to-back paper [122] (with extended version [123]), the results for the Néel quench were reproduced and extended to the dimer product state

$$\begin{eqnarray}&&|\Psi(t=0)\rangle =\frac{1}{{{2}^{(N+1)/2}}}(|(\uparrow \downarrow -\downarrow \uparrow )(\uparrow \downarrow -\downarrow \uparrow )...\rangle \,+\,|(\downarrow \uparrow -\uparrow \downarrow )(\downarrow \uparrow -\uparrow \downarrow )...\rangle )\end{eqnarray} \tag{ 129 }$$

(corresponding to the ground state of the Majumdar–Ghosh model [124]). The great virtue of the dimer state is that it showed an even more dramatic contrast between the Quench Action and GGE predictions than the Néel case, with the discrepancies being clearly exposed in numerical verifications using infinite size time evolved block decimation. This study was later extended to include a wider class of q-dimer initial states (for which exact overlaps are known [112]) in [123], including more a detailed study of the local longitudinal and transverse spin–spin correlations. These studies provided additional confirmation that the Quench Action provided the correct answers for the steady state.

As far as the steady state is concerned, the story concerning the discrepancy between the exact Quench Action results and the GGE was brought to a close with the understanding that the correct GGE [125] in the case of quantum spin chains needed to also involve families of quasilocal conserved charges generated from higher-spin transfer matrices. We refer to the accompanying paper [126] for an extensive introduction to these quasilocal charges and to their use in different contexts, and to [127] for a rigorous discussion of locality as far as thermalization is concerned. We however point out that these quasilocal charges, first derived in [128], can be used when quenching to the gapless regime [129], this representing another road to obtain the steady state (but not the time evolution).

To close this section, although it is not our intent to review the field here, let us nonetheless briefly mention that the Quench Action has already been used in many other contexts in addition to those already described, including the sine-Gordon [130] and sinh-Gordon [131] models, geometric quenches in free fermionic chains [132], finite integrable spin chains [133], the Kondo model [134], the Lieb–Liniger gas with hard-wall boundary conditions [135], Bragg pulses in bosonic gases [136], quenches from rotating BECs [137], etc. Another very interesting case is that of the extension of the BEC to repulsive Lieb–Liniger quench work [36] to the attractive case [138], which also relates to the KPZ equation [139]. Since this model supports bound state excitations, the implementation of the Quench Action requires treating multiple root distributions, similarly to what needs to be done in the case of quantum spin chains. These illustrate the variety of contexts in which the method can be fruitfully applied.

Taking a step away from specific cases of implementation of the Quench Action, it is worthwhile discussing the most important building blocks of the Quench Action, namely the exact overlaps which stand at the base of the whole formalism. Restricting our discussion for the moment to the case of quantum quenches, the general problem is to find workable expressions for the overlaps

$$\begin{eqnarray}&&{{\langle {{\Psi}_{0}}|\left\{I\right\}\rangle}_{g}}\equiv {{\langle {{\Psi}_{0}}|\left\{\lambda \right\}\rangle}_{g}}\end{eqnarray} \tag{ 130 }$$

between some initial state $|\Psi\rangle$ and the eigenstates of a certain (preferably, though not necessarily integrable) Hamiltonian with interaction parameter g, which we have denoted alternately in terms of quantum numbers or rapidities solving the relevant Bethe equations (including normalization).

The two cases we have discussed in more detail, namely the BEC to Lieb–Liniger quench on the one hand, and the Néel to XXZ quench in the other, rely on overlaps described in [110, 111, 113, 114]. Other known results include overlaps between states at $\Delta =\infty$ and $\Delta =0$ [140] and of partial Néel states [141]. An interesting further proposal provides some recursive expressions for overlaps of simple product states [142]. Perhaps one could extend these all the way to any initial state satisfying cluster decomposition, these being representative of the generic case (see for example ([143]).

It is however not trivial to generalize these results to more general initial states. In fact, these two problems share many similarities with each other. The resemblance between (95) and (121) is particularly striking, and we can wonder how general similar constructions could be. Gaudin's formula, giving the norm of a Bethe eigenstate in terms of the determinant of the Gaudin matrix, has a clear intuitive origin: the quantum numbers of an eigenstate can be viewed as the proper independent variables labelling a state, and the measure for summation over eigenstates is flat in quantum number space. Since the Bethe equations define the mapping between quantum number space and rapidity space, the measure for summation over eigenstates, when translated to rapidity space, gets rescaled by the Jacobian of the transformation between quantum numbers and rapidities, namely by the determinant of the Gaudin matrix whose entries are the derivatives of the Bethe equations (second derivatives of the Yang–Yang action). The fact that the overlaps (95) and (121) are given by a Gaudin-like matrix is however less intuitive from the outset: these overlaps can be viewed as a kind of 'square-root norm' originating from the partition function of the 6-vertex model with reflecting boundary conditions [110], but the deeper meaning of these identities remains obscure at this stage.

An interesting open challenge is to try to generalize (95) to the case of a generic free boson initial state

$$\begin{eqnarray}&&|{{\Psi}_{0}}\rangle =\,|{{\left\{k\right\}}_{N}}\rangle.\end{eqnarray} \tag{ 131 }$$

For generic choices of momenta {k}_N, the overlap with Bethe states cannot be of Gaudin-like form (95). In fact, the perspective for the obtention of a genuinely useful representation of the overlaps is pessimistic: one can create the states (131) by applying a product of single-particle creation operators on the vacuum; alternately, the overlap can be viewed as the overlap of the vacuum with the remnants of the Bethe states after acting with the free boson annihilation operators. Whether the in principle factorially large sum which remains can be expressed in an economical determinant-like form remains an open issue at this stage. Note however that the obtention of a workable formula here would in principle solve the extremely interesting and important problem of calculating overlaps of eigenstates at different values of the interaction parameter,

$$\begin{eqnarray}&&{{}_{g}}{{\langle \left\{I\right\}|\left\{I\right\}\rangle}_{{{g}^{\prime}}}},\end{eqnarray} \tag{ 132 }$$

by simply using the free bosonic basis as an intermediate basis. We leave this problem as an outstanding challenge to the community, noting in passing that there is one case in which such a formula can be obtained: for the Richardson model

$$\begin{eqnarray}&&H=\underset{\alpha}{\sum}\underset{\sigma =\uparrow,\downarrow}{\sum}\frac{{{\varepsilon}_{\alpha}}}{2}c_{\alpha \sigma}^{\dagger}{{c}_{\alpha \sigma}}-g\underset{{{\alpha}_{1}},{{\alpha}_{2}}}{\sum}c_{{{\alpha}_{1}}\uparrow}^{\dagger}c_{{{\alpha}_{1}}\downarrow}^{\dagger}{{c}_{{{\alpha}_{2}}\downarrow}}{{c}_{{{\alpha}_{2}}\uparrow}},\end{eqnarray} \tag{ 133 }$$

one can write the exact overlap between states at different values of g [144] by invoking Slavnov's theorem. Although the Quench Action method is not usable here (the thermodynamic limit of this model is inherently mean-field), the algebraic structure of the overlaps might serve as inspiration for other non mean-field-like cases.

The problem we just mentioned, namely that of the overlap between exact eigenstates of Hamiltonians with differing interaction parameters, is formally written in terms of algebraic Bethe Ansatz operators as

$$\begin{eqnarray}&&\langle 0|\underset{j}{\prod}{{C}^{(g)}}(\lambda _{j}^{(g)})\underset{k}{\prod}{{B}^{({{g}^{\prime}})}}(\mu _{k}^{({{g}^{\prime}})})|0\rangle \end{eqnarray} \tag{ 134 }$$

in which we have explicitly labeled the monodromy matrix operators B, C in terms of the interaction parameter since the algebra of these operators explicitly depends on these interaction parameters. Similarly, the sets of rapidities are superscripted by the value of the interaction parameter for which they solve Bethe equations. Although the commutation relations between ABA operators pertaining to a given interaction parameter are relatively simple (namely the usual quadratic algebra), the 'cross'-algebra giving the commutation relations between operators at different interaction values is dramatically more complicated. One can however dream of a higher algebraic structure intertwining the ABA operators at g with those at ${{g}^{\prime}}$; this structure, which might take the form of some exponentiated, continuous unitary transformation, remains however elusive at this stage. A more immediately workable but still interesting possibility is thus to exploit the crucial feature making Slavnov's theorem interesting, namely that the overlaps

$$\begin{eqnarray}&&\langle 0|\underset{j}{\prod}C({{\lambda}_{j}})\underset{k}{\prod}B({{\mu}_{k}})|0\rangle \end{eqnarray} \tag{ 135 }$$

(in which we now consider a single interaction parameter g) are exactly known provided one of the sets $\{\lambda \}$ or $\{\mu \}$ solves the Bethe equations, the other set being completely arbitrary. It is thus possible to immediately write ovelaps between Bethe states and generic initial states

$$\begin{eqnarray}&&|{{\Psi}_{0}}\rangle \equiv |\{\mu \}\rangle \end{eqnarray} \tag{ 136 }$$

for arbitrary (but judiciously chosen so at to make the problem interesting) state-defining rapidity set $\left\{\mu \right\}$. Although the precise nature of such states is not immediately obvious, it can in principle be quantified by mapping the states back to Bethe states (again invoking Slavnov's theorem) and computing representative observables. By cleverly choosing these sets of rapidities based on known intuitions, one can define whole families of quench problems whose solution via the Quench Action becomes straightforward. We suggest this as a line worth exploring in the future.

Finally, an important point to bear in mind for the application of the Quench Action method is that the knowledge of the exact overlaps to all orders in system size is overkill, at least as far as the determination of the steady state is concerned. For this purpose, only the leading extensive part of the logarithmic overlap is needed. It might be possible to devise simplified schemes in which such leading parts can be obtained, bypassing exact calculations. Results on these, when fed back into the Quench Action formalism, would allow to reconstruct the full time evolution of observables from quench time to the steady state. We again leave this as an open challenge.

In both the BEC to Lieb–Liniger and the Néel to XXZ quenches, remarkably, the simplicity of the original state has imprinted itself in the analyticity of the solution to the GTBA for the steady state. The existence of a closed-form solutions for (generalized) TBA equations is atypical (an interesting early example being the appearance of Airy functions in integrable N  =  2 supersymmetric models [145]). What the deeper meaning of this fact is, and whether such a 'simplicity-analyticity' correspondence holds for other quenches, are interesting questions for the future.

The one direction which probably offers the most prospects for new results is the one of the analysis of the time-dependent behaviour of observables after quenches. It is natural for theoretical efforts to have been up to now mostly directed at the characterization of steady states, since off-diagonal terms in the expectation values can then be dropped. On the other hand, the presence of any integrability-breaking term [146–150], even arbitrarily weak, destroys the long-time limit, presumably thermalizing it, though the timescales for this can be very large. Shifting the focus away from the steady state, it is obvious that even more interesting physics is to be obtained in the actual relaxation process itself, for which the Quench Action offers a handle. We have showed some results on this in the case of the BEC quench in Lieb–Liniger. These constructions could be extended to more complicated observables at the cost of surmountable computational difficulties. In the Tonks–Girardeau limit, the explicit time dependence can be constructed (from direct calculations which are easily reproducible from the Quench Action, see [136]) for quenches such as the trap release [51, 52] (also to a gas in a hard-wall trap [151]) or a Bragg pulse [136]. In the low-density limit, generic quenches also lead to Tonks–Girardeau behaviour [49].

Time dependence in the case of the Néel quench or of the dimer quench in spin chains (studied in e.g. [119]) is something which is in principle accessible to the Quench Action but remains to be done. At generic interaction value, for complicated but meaningful quantities such as (relative) time- and space-dependent spin correlations, one difficulty is that it then requires the computation of matrix elements on states with mixed string contents, in highly-excited quantum number configurations. This has been explored in some recent works [63, 80, 152]. Numerically, it is possible to perform the necessary summations in (66) by using the ABACUS algorithm, whose logic is applicable starting from any (even highly excited) state, though the summations become difficult to perform due to the finite entropy of the starting state and to the large system sizes required.

We have seen before that the implementation of the Quench Action requires as ingredients: the saddle-point distribution ${{\rho}_{\text{sp}}}$ and the characteristic quench overlap function s. Given the excitation dispersion relations ε (defined by ${{\rho}_{\text{sp}}}$ if the Hamiltonian is specified), and given an observable $\mathcal{O}$ which we are interested in, a knowledge of its matrix elements then allows to compute the post-quench time dependent expectation value, using the fundamental representation (66), 'working backwards' from the steady state by including more and more excitations.

Going further with the idea of 'working backwards', one can simply turn the tables around and back-engineer quench problems by postulating the saddle-point ${{\rho}_{\text{sp}}}$ and characteristic overlap function s. Given these, one can then try to solve the inverse problem of determining which kind of state stood at the pre-quench origin of time. The level of exactitude with which this state can be determined depends on the detailed knowledge of s. Starting with the basics, one can imagine that many different initial states relax to the same saddle-point ${{\rho}_{\text{sp}}}$. How this relaxation occurs then becomes dependent on the interplay between the characteristics of ${{\rho}_{\text{sp}}}$ (e.g. its inflection points) and those of the overlap function s and dispersion relation ε. One could thus develop the notion of 'universality classes' of types of relaxation behaviour based on these characteristics.

An interesting exploratory route is to consider a much more brutally out-of-equilibrium situation, namely a periodically-driven system, in other words Floquet dynamics. The driving would then be periodic over a period T; in its simplest form, one can envision a 'quench-dequench' cycle switching between Hamiltonians H₁ (for a duration of t₁) and H₂ (duration t₂, with $T={{t}_{1}}+{{t}_{2}}$) with Floquet operator

$$\begin{eqnarray}&&{{\text{e}}^{-\text{i}{{t}_{2}}{{H}_{2}}}}\,{{\text{e}}^{-\text{i}{{t}_{1}}{{H}_{1}}}}.\end{eqnarray} \tag{ 137 }$$

Starting from a given initial state, repeated application of the Floquet operator will keep the system from equilibrating. If one focuses on observables measured at stroboscopic times (namely always at the same instant in the Floquet cycle), one can however expect pseudo-relaxation, in the sense that the system would tend to effectively relax to one of the set of Floquet eigenstates. The logic of the Quench Action can be applied to describe such stroboscopic observables, the Floquet eigenstates now taking on the previous role played by the Hamiltonian eigenstates in the quench problem. A steady-state can be expected (i.e. relaxation to a saddle-point Floquet eigenstate, see an example in spin chains in [153]) again due to the large relative dephasings. Finding this stroboscopic steady state would again be the result of an optimization procedure following the logic of the Quench Action. We have not yet been able to control such a calculation, and leave the finding of such an example as another open problem.

Although we have used the language of Bethe Ansatz to present the logic of the Quench Action, this is not a necessary requirement. Given a computationally useful basis of eigenstates, precisely the same steps can be followed in complete generality, the whole edifice being based on the sole existence of such an eigenbasis. The reason why Bethe Ansatz solvable models are productively used in actual implementations of the Quench Action is that they provide all the necessary ingredients: states and their quantum number labels, overlaps, matrix elements. One could imagine for example a future in which numerical methods, perhaps based on matrix-product states, are able to deliver a set of practical tools from which quench problems can be treated using the Quench Action logic. Admittedly, this is at the moment wishful thinking, but who knows what the future will bring.