Out-of-Equilibrium Dynamics of Strongly-Correlated Systems:
Lessons from Integrability

Séminaire LPMMC

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

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 the full time evolution of observables as

$ \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