The ins and outs of quantum dynamics
Integrability,
from Newton's cradle to Gibbs' grave

Polish Academy of Sciences

Jean-Sébastien Caux

Institute of Physics

University of Amsterdam

Out-of-equilibrium: why?

Why NOT?

• Real systems equilibrate very rapidly
• Experimentally: tough to control
• Theoretically: can't do it - toolbox too limited

Why INDEED?

• Real isolated systems can equilibrate slowly
• Novel states & properties
• Experimentally: frontier of current capabilities
• Theoretically: the toughest challenge around?

Out of equilibrium: categories

• local, gentle (linear response)
• local, hard (nonlinear response)
• local, gentle, repeated
• local, hard, repeated
• local, gentle, erratic
• local, hard, erratic
• extended, one-time
• extended, repeated regularly
• extended, erratic

Out of equilibrium: terminology

Out of equilibrium: what? how?

Out-of-equilibrium using Integrability

The Lieb-Liniger Model

$ H = \int_0^L dx \left[ \partial_x \Psi^{\dagger}(x) \partial_x \Psi(x) + c \Psi^{\dagger} (x) \Psi^{\dagger}(x) \Psi(x) \Psi(x) \right] $

The Lieb-Liniger Model

The Lieb-Liniger Model

Description using the Algebraic Bethe Ansatz

\[ \tau(\lambda) = tr T(\lambda), \hspace{5mm} T(\lambda) = \left( \begin{array}{cc} A(\lambda) & B(\lambda) \\ C(\lambda) & D(\lambda) \end{array} \right), \hspace{5mm} | \Psi_N \left( \left\{ \lambda \right\}_N \right \rangle = \prod_j B(\lambda_j) | 0 \rangle \]

Matrix elements of physical operators known from:

• Gaudin formula for state norm

  \[ \langle \{ \lambda \} | \{ \lambda \} \rangle = c^N \prod_{j > k} \frac{\lambda_{jk}^2 + c^2}{\lambda_{jk}^2} \mbox{det} \left[ \delta_{jk} \left(L + \sum_l K_{jl} \right) - K_{jk} \right], \hspace{10mm} K_{jk} = \frac{2c}{\lambda_{jk}^2 + c^2} \]

• solution of inverse problem $A, B, C, D \leftrightarrow \psi(x), \psi^\dagger(x)$
• Slavnov formula for $\langle 0 | \prod_j C(\lambda_j) \prod_k B(\mu_k) | 0 \rangle$

The Lieb-Liniger Model

(Dynamical) Correlation functions

\[ {\cal S}^{a \bar{a}} (k, \omega) \equiv \frac{1}{N} \sum_{j, j'} e^{-i k(j-j')} \int_{-\infty}^\infty dt e^{i\omega t} \langle \frac{1}{2} \left[ {\cal O}^a_j (t), ({\cal O}^a_{j'} (0))^\dagger \right] \rangle \]

Simplest examples:

• Dynamical structure factor: ${\cal O}_j = \rho_j \equiv \Psi^\dagger(x_j) \Psi(x_j)$
• One-body function (Green's): ${\cal O}_j = \Psi(x_j)$ or $\Psi^\dagger(x_j)$

Dynamical structure factor
from (algebraic) Bethe Ansatz

Abacus (post Achilles)

Dynamical structure factor
from (algebraic) Bethe Ansatz

Lieb-Liniger using cold atoms (I)

Lieb-Liniger using cold atoms (II)

Interaction quench in Lieb-Liniger

We will consider a sudden switch of the interaction:

$ c (t) = \left\{ \begin{array}{ll} c_i & t < 0, \\ c_f & t > 0 \end{array} \right. $

starting in ground/(any eigen)state of $H(c_i)$ at $t < 0$

For many reasons, this is a complicated problem...

• GGE logic fails, charges have $\infty$ exp. values
• approaches using regularizations (q-bosons, extended charges, ...): no complete solution

The strategy

We write the initial state in the basis of the eigenstates of $H(c_f)$, denoted as $\left\{ |\{\lambda\}^{(n)}_N\rangle \right\}, n = 0, 1, ...$

Setup

The initial state is decomposed in the basis of post-quench Hamiltonian eigenstates:

$ |\Psi (t=0) \rangle = \sum_{\{ I \}} e^{-S^\Psi_{\{ I \}}} | \{ I \} \rangle, $

$ S^\Psi_{\{ I \}} = -\ln \langle \{ I \} | \Psi (t=0) \rangle ~\in ~\mathbb{C} $

Going to the thermodynamic limit

Time-dependent expectation values

In ThLim, saddle-point evaluation $\left. \frac{\delta S^Q[\rho]}{\delta \rho} \right|_{\rho_{sp}} = 0$ exactly gives full time evolution of observables:

$ \begin{align} \lim_{Th,reg} \bar{\cal O} (t) = \lim_{Th,reg} \frac{1}{2} \sum_{\bf e} \left[ e^{-\delta S_{\bf e}[\rho_{sp}] - i \omega_{\bf e}[\rho_{sp}] t} \langle \rho_{sp} | {\cal O} | \rho_{sp}; {\bf e} \rangle \right. \nonumber \\ \left. + e^{-\delta S^*_{\bf e}[\rho_{sp}] + i \omega_{\bf e}[\rho_{sp}] t} \langle \rho_{sp}; {\bf e} | {\cal O} | \rho_{sp} \rangle \right] \end{align} $

The BEC to Lieb-Liniger quench

Initial state: ground state of free bosons

$ | \Psi_0 \rangle = \frac{1}{\sqrt{L^N N!}} \left(\psi_{k=0}^\dagger\right)^N | 0 \rangle $

Exact solution using Quench Action

$ S^Q = L\int_0^\infty d\lambda \rho(\lambda) \log\left(\frac{\lambda^2}{c^2} \left( \frac{1}{4} + \frac{ \lambda^2}{c^2} \right) \right) + S^{YY} $

Exact solution using Quench Action

Steady state distribution: non-thermal
J. De Nardis, B. Wouters, M. Brockmann & J-SC, PRA 89, 2014

Viewed as a perturbation

$ H(c_i) = H(c_f) + (c_i - c_f) \int_0^L dx \Psi^\dagger(x)\Psi^\dagger(x)\Psi(x)\Psi(x) $

Matrix elements of the original Hamiltonian in the computational basis of post-quench eigenstates:

$ \begin{align} \langle \{\lambda\}^{(m)}_N | H(c_i) | \{\lambda\}^{(n)}_N \rangle &= \delta_{n,m} E\big(\{\lambda\}^{(n)}_N \big) + \\ &+ (c_i - c_f) L\, \langle \{\lambda\}^{(m)}_N | \big(\Psi^\dagger(0)\big)^2 \big( \Psi(0) \big)^2 |\{\lambda\}^{(n)}_N\rangle \end{align} $

Said otherwise, we perturb the system
with the operator

$ g_2(0) = \big( \Psi^\dagger(0)\big)^2 \big(\Psi(0)\big)^2 $

Truncated Spectrum Approach

Idea: truncate Hilbert space to low-E sector using cutoff $\Lambda$, then diagonalize numerically

NRG + TSA

The overall workflow:

• Construct the computational basis $\{ |\{\lambda\}^{(j)}\rangle\}$ and order by energy $E(\{\lambda\}^{(j)})$
• Construct the Hamiltonian with the lowest $N_s+\Delta N_s$ states: $ \begin{align} \Big\{ |\{\lambda\}^{(1)}\rangle, \ldots, |\{\lambda\}^{(N_s+\Delta N_s)}\rangle\Big\}, \end{align} $ and diagonalize it, obtaining approximate energies and eigenstates $ \begin{align} \Big\{ |E^{(1)}\rangle, \ldots , |E^{(N_s+\Delta N_s)}\rangle \Big\}. \end{align} $
• Discard the highest $\Delta N_s$ approximate eigenstates $ \begin{align} \Big\{|E^{(N_s+1)}\rangle,\ldots,|E^{(N_s+\Delta N_s)}\rangle\Big\}, \end{align} $ and form a new basis with the lowest $N_s$ approximate eigenstates and the next $\Delta N_s$ lowest states in the computational basis.
• Construct the Hamiltonian in the new basis, and diagonalize to obtain new approximate eigenstates.
• Return to the third step.

Harmonically trapped Lieb-Liniger

JSC & R. M. Konik, Phys. Rev. Lett. 109, 175301 (2012)

Interaction quench: GSE from NRG + TSA

Neil Robinson, Bart de Klerk + JSC,  SciPost Phys. 11, 104 (2021)

Ground-state energy from NRG + TSA

NRG follows TSA, but can include many more states.

Is ordering in energy a good idea?

Let's try an alternative metric...

We will order by matrix element size, between states in computational basis and the 3 lowest-lying ones from TSA:

$ \left| \langle \{\lambda\}^{(n)}_N | g_2(0) | \tilde E_j\rangle \right|, ~~j =0,1,2. $

NRG-E versus NRG-ME

Energy ordering versus Matrix Element ordering

NRG-ME: convergence of the overlaps

Towards a "perfect" ordering

ABACUS is really efficient at computing
dynamical correlation functions

For an operator ${\cal O}$, it preferentially seeks states with high matrix elements $ \langle \{\lambda\}^{(0)}_N | {\cal O} | \{\lambda\}^{(n)}_N\rangle $

HOSTS: High Overlap States Truncation Scheme

This produces an ordering which is similar to the "back-engineered" one from the overlaps:

Energy convergence using HOSTS

• ABACUS short run
• 390K states generated
• preferentially ordered
• NRG performed with $N_{NRG}=800$ and $\Delta N = 160$
• Convergence of $E_0$ to under 1% achieved with only a few thousand basis states

Convergence of other quantities using HOSTS

Overlaps

Local expectation value: $g_2$

Real-time dynamics after the quench

Return amplitude

$ \langle \Psi_i |\Psi_i(t)\rangle \approx \sum_{n=0}^{N_{\text{tot}}} e^{-\mathrm{i}E(\{\lambda\}^{(n)}t} \left\vert \langle \{\lambda\}^{(n)} | \Psi_i \rangle \right\vert^2 $

Fidelity

$ {\cal F}(t) = |\langle\Psi_i|\Psi_i(t)\rangle|^2 $

Time dependence of local observables

$ \begin{align} \langle &O(t)\rangle_i \equiv \langle \Psi_i(t) | O |\Psi_i(t)\rangle\\ & = \sum_{n,m=0}^{N_\text{tot}} e^{-\mathrm{i}t[E(\{\lambda\}^{(n)}) - E(\{\lambda\}^{(m)})]} \langle \Psi_i | \{ \lambda \}^{(m)} \rangle \langle \{ \lambda \}^{(m)} | O |\{\lambda\}^{(n)} \rangle \langle \{\lambda\}^{(n)} | \Psi_i\rangle \end{align} $

Although pretty successful, our alternative ordering still somehow relies on the perturbation being weak.

Desire: describe strongly non-perturbative quenches.

Challenge: account for contributions from important intermediate states.

Idea: keep approximate states at each step of the NRG, and select at each individual step which ones should participate, based on quadratic approximation of perturbative series.

$\rightarrow$ Matrix Element Renormalization

Matrix Element Renormalization

The overall workflow (steps 1-3 / 7):

• Generate the computational basis via preferential state generation from the seed state $|\Omega \rangle$. Order the states in the computational basis according to the metric "ME over Ediff" to obtain $\{|\{\lambda\}^{(j)}\rangle \}$.
• Construct a truncated Hamiltonian from the first $N_s+\Delta N_s$ computational basis states, $\big\{ |\{\lambda\}^{(1)}\rangle, \ldots, |\{\lambda\}^{(N_s+\Delta N_s)}\rangle\big\}$ and diagonalize this Hamiltonian to obtain the first approximate eigenstates $\big\{|1\rangle, \ldots , |N_s+\Delta N_s\rangle\big\}$ with energies $\{ E_1, \ldots, E_{N_s + \Delta N_s} \}$.
• Compute the overlaps between $|\{\lambda\}^{(0)}\rangle$ and the newly acquired approximate eigenvectors $\{|1\rangle, \ldots, |N_s + \Delta N_s\rangle \}$. Then relabel the approximate eigenstates such that $|1\rangle$ refers to the approximate eigenstate with the largest overlap.

Matrix Element Renormalization

The overall workflow (step 4 / 7):

• Take the next $\Delta N_s$ computational basis states $|\{\lambda\}^{(j)}\rangle$ and compute the second order weight for each of the approximate eigenstates $|i \rangle$ with $i>1$ obtained in previous steps: \[ w_2\left( |i \rangle \right) = \sum_{j} \frac{ \langle i | \delta H |\{\lambda\}^{(j)} \rangle \langle \{\lambda\}^{(j)} | \delta H | 1 \rangle}{( E_1 - E_{\{\lambda\}}^{(j)}) (E_1 - E_i)}. \] Here $\delta H = cR\times g_2(0)$ is the perturbing operator, the sum ranges over all the newly added computational basis states, and $E_{\{\lambda\}}^{(j)} = E(\{\lambda\}^{(j)})$ are the energies of the newly added computational basis states.

Matrix Element Renormalization

The overall workflow (steps 5 - 7 / 7):

• Form a truncated basis consisting of $|1\rangle$, the $N_s-1$ states in $\{|i \rangle \big\}_{i>1}$ with the largest $w_2$-weight, and the $\Delta N_s$ computational basis states introduced in step 4, construct the truncated Hamiltonian in this basis and diagonalize it to obtain $N_s + \Delta N_s$ new approximate eigenstates $\{|1'\rangle, \ldots, |(N_s + \Delta N_s)'\rangle\}$.
• Compute the overlaps between $|1\rangle$ and the newly acquired approximate eigenvectors, and replace $|1\rangle$ by the approximate eigenvector with the largest overlap. Replace the remaining approximate eigenvectors used to form the truncated basis in step 5 with the remaining newly obtained approximate eigenvectors.
• Return to the fourth step.

Convergence using MERG

Overlaps

Energy

MERG for a harder quench

Quenches from spatially inhomogeneous states

After initial dephasings: 'hydrodynamic' evolution described by local GGE

$\partial_t \rho (\lambda) + \partial_x (v^{\tiny \mbox{eff}}(\lambda) \rho (\lambda)) = 0$

with effective velocities determined by the dressing equation

$v^{\tiny \mbox{eff}}(\lambda) = v^{\tiny \mbox{gr}}(\lambda) + \int d\lambda' \frac{\varphi(\lambda, \lambda')}{p'(\lambda)} \rho(\lambda') \left(v^{\tiny \mbox{eff}}(\lambda') - v^{\tiny \mbox{eff}}(\lambda)\right)$

Effective dynamics:
"flea gas" of scattering quasiparticles

Flea gas GHD
for the Quantum Newton's Cradle

JSC, B. Doyon, J. Dubail, R. Konik and T. Yoshimura, SciPost Phys. 6, 070 (2019)

Short time scales:

Flea gas GHD
for the Quantum Newton's Cradle

JSC, B. Doyon, J. Dubail, R. Konik and T. Yoshimura, SciPost Phys. 6, 070 (2019)

Longer time scales:

Generalized Hydrodynamics with space/time-dependent potentials

Phys. Rev. Lett. 123, 130602 (2019)

Adiabatic formation of bound states in the 1d Bose gas

Phys. Rev. B 103, 165121 (2021)

Observing Bethe strings in an attractive Bose gas far from equilibrium

Milena Horvath, Alvise Bastianello, Sudipta Dhar, Rebekka Koch, Yanliang Guo,
Jean-Sébastien Caux, Manuele Landini, and Hanns-Christoph Nägerl, arXiv:2505.10550

Observing Bethe strings in an attractive Bose gas far from equilibrium

Milena Horvath, Alvise Bastianello, Sudipta Dhar, Rebekka Koch, Yanliang Guo,
Jean-Sébastien Caux, Manuele Landini, and Hanns-Christoph Nägerl, arXiv:2505.10550

The structure factor
of the Heisenberg chain

Quantum wake dynamics in Heisenberg

A. Scheie, P. Laurell, B. Lake, S. E. Nagler, M. B. Stone, JSC, D. A. Tennant
Nature Communications 13, 5796 (2022)

• Inelastic neutron scattering on $K Cu F_3$
• Really accurate data in momentum/frequency space enabling...
• Fourier transform to real space / time
• Look specifically at high-energy features

Quantum wake dynamics in Heisenberg

A. Scheie, P. Laurell, B. Lake, S. E. Nagler, M. B. Stone, JSC, D. A. Tennant
Nature Communications 13, 5796 (2022)

• Quantum wake:
  a local, coherent, long-lived, quasiperiodically oscillating magnetic state emerging out of the distillation of propagating excitations following a local quantum quench

Quantum wake dynamics in Heisenberg

Quantum wake dynamics in Heisenberg

Quantum wake dynamics in Heisenberg

Selectively (patch-wise) Fourier transforming: $k=\pi/2$ (mod $2\pi$) states are responsible for the quantum wake dynamics

Beyond the reach of any low-energy method