Opening Krylov space to access all-time dynamics via dynamical symmetries

Basic Information

• Paper ID: 2503.07403
• Title: Opening Krylov space to access all-time dynamics via dynamical symmetries
• Authors: Nicolas Loizeau, Berislav Buča, Dries Sels
• Institutions: Niels Bohr Institute (Copenhagen), Université Paris-Saclay, University of Oxford, New York University, Flatiron Institute
• Classification: quant-ph (Quantum Physics)
• Publication Date: August 22, 2025 (v2 version)
• Paper Link: https://arxiv.org/abs/2503.07403v2

Abstract

This paper addresses the fundamental challenge of understanding short-time and long-time dynamics in closed quantum many-body systems. For closed systems with local interactions, the dynamics of local observables can be expanded in terms of (pseudo-local) eigenmodes of the Liouvillian, known as dynamical symmetries. These operators fall into two categories: transient operators that decay with time, and eternal operators that oscillate permanently or remain invariant (conservation laws). The paper proposes a method to numerically and analytically derive these dynamical symmetries by introducing naturally emergent open boundary conditions on the Krylov chain, partitioning the Krylov space into system and environment degrees of freedom, where non-local operators constitute an effective bath for local operators.

Research Background and Motivation

Core Problems

1. Computational Challenges in Quantum Many-Body Dynamics: Solving short-time and long-time dynamics of closed quantum many-body systems represents one of the primary challenges in atomic and condensed matter physics
2. Difficulty in Extracting Dynamical Symmetries: Although theoretically the dynamics of local observables can be completely characterized by eigenmodes of the Liouvillian (dynamical symmetries), deriving these operators has been impossible except for a few special models

Problem Significance

• Complete Dynamical Characterization: Dynamical symmetries provide complete characterization of system dynamics. For any operator O and initial state ρ, the time evolution can be expressed as: Tr(O(t)ρ) = ∑_ω e^(iωt)μ_ω
• Understanding Non-equilibrium Phenomena: Essential for understanding many-body localization, time crystals, quantum scars, fragmentation, and other non-trivial dynamical behaviors
• Thermalization Mechanisms: In chaotic systems, connects the growth of correlations to thermalization and exponential decay of observables

Limitations of Existing Methods

1. Finite System Constraints: The Liouvillian spectrum of finite closed systems is purely real, preventing access to thermodynamic limit properties
2. Pseudo-mode Expansion Methods: Previous approaches (Refs 96,97) are limited to dynamics near equilibrium
3. Pseudo-locality Unexplored: Prior studies did not examine the (pseudo-)locality of dynamical symmetries, which is crucial for non-equilibrium physics

Research Motivation

Propose a more rigorous method that achieves system-environment decomposition by introducing open boundary conditions on the Krylov chain, enabling:

• Local quantities supported in the system
• Non-local quantities forming the environment
• Identification of which dynamical symmetries are pseudo-local and therefore relevant

Core Contributions

1. Novel Open Boundary Condition Method: Proposes introducing dissipative boundary conditions at a specific position L on the Krylov chain, avoiding non-physical wave function reflections, with the form: ∂t φ_L = (b_L + b{L+1})φ_ - 2b_{L+1}φ_L
2. System-Environment Decomposition Framework: Establishes a natural decomposition of Krylov space where:
   • Left side (system): simple local operators with exactly known hopping coefficients
   • Right side (environment): complex non-local operators forming an effective bath
3. Analytical Solutions for Ideal Cases: Provides exact solutions for two idealized scenarios:
   • Linear growth b_n = n: eigenvectors expressed in terms of Meixner polynomials
   • Square-root growth b_n = √n: eigenvectors expressed in terms of Hermite polynomials
4. Direct Connection Between Operator Growth and Thermalization: Proves that in the linear growth case (chaotic systems), all observables decay exponentially at rate -2α, where α is the growth rate of Lanczos coefficients
5. Numerical Verification and Applications: Validates the method on actual spin systems such as XXZ chains and transverse-field Ising chains, successfully extracting dynamical symmetries

Detailed Methodology

Task Definition

Input:

• Hamiltonian H
• Initial operator O_0 (or initial state density matrix ρ)

Output:

• Dynamical symmetries A_ω satisfying H, A_ω = -ωA_ω
• Frequencies ω (can be complex): real part corresponds to oscillations, imaginary part to decay

Constraints:

• Focus only on pseudo-local dynamical symmetries (affecting local physics)
• Applicable to systems with local interactions and short-range correlated initial states

Model Architecture

1. Krylov Space Construction (Lanczos Algorithm)

Starting from seed operator O_0, recursively apply the Liouvillian L = H, · to construct an orthonormal basis {O_n}:

O_1 = LO_0/b_1 = [H,O_0]/b_1,  b_1 = ||LO_0||

For n > 2:
O'_n = LO_{n-1} - b_{n-1}O_{n-2}
O_n = O'_n/b_n
b_n = ||O'_n||

where the norm is defined as ||O||² = 1/(2N) TrO²

2. One-Dimensional Single-Particle Mapping

Time evolution in the Krylov basis maps to a 1D single-particle hopping problem:

O(t) = 1/(2N) ∑_n i^n φ_n(t)O_n

Coefficient evolution equation: ∂t φ_n = b_n φ - b_{n+1} φ_{n+1}, φ_n(0) = δ_

3. Open Boundary Conditions

Truncate the chain at position L, assuming φ_n is sufficiently smooth to be locally approximated as linear: φ_{L+1} ≈ φ_L + (φ_L - φ_)

Yielding the dissipative boundary condition: ∂t φ_L = (b_L + b{L+1})φ_ - 2b_{L+1}φ_L

4. Modified Liouvillian Matrix

The truncated Liouvillian has a tridiagonal plus boundary correction form:

L̃ = i × [tridiagonal matrix, last row special form]
Last row: [0, 0, ..., b_L + b_{L+1}, -2b_{L+1}]

This non-Hermitian matrix allows complex eigenvalues

Technical Innovations

1. Physical Meaning of Boundary Conditions

• No Parameter Introduction: Unlike Refs 96,97 with external dissipation + extrapolation, this method requires no additional parameters
• Smoothness Assumption: Based on natural assumption of derivative continuity of Krylov wave functions
• Local-Nonlocal Separation: Krylov chain right-side operators are dominated by k>L k-local Pauli strings

2. Distinction from Traditional Methods

• Dirichlet Boundary: φ_{L+1} = 0 leads to non-physical reflection
• This Method: Allows flow into environment and return, more realistically reflecting operator propagation

3. Design Rationality

In ideal cases (b_n = n or √n):

• Dynamics under open boundary conditions nearly indistinguishable from infinite chain (Figure 2)
• Real part of spectrum nearly unchanged, imaginary part introduces physical decay (Figure 3)
• Linear growth case: all non-trivial roots satisfy Im(ω) = -2i (numerically verified to n=2^12)

Experimental Setup

Model Systems

1. XXZ Chain (Integrable Model)

H = ∑_i (s^x_i s^x_{i+1} + s^y_i s^y_{i+1} + Δs^z_i s^z_{i+1} + hs^z_i)
Parameters: Δ = -1/2, h = 2
Initial operator: O_0 = Q_3 = ∑_i (s^+_i s^+_{i+1}s^+_{i+2} + s^-_i s^-_{i+1}s^-_{i+2})

Known to possess dynamical symmetry with frequency ω = 12

2. Chaotic Chain

H = ∑_i (s^x_i s^x_{i+1} - 1.05s^z_i + 0.5s^x_i)
Initial operator: O_0 = ∑_i (1.05s^x_i s^x_{i+1} + s^z_i)

Operator chosen with no overlap with H to exclude trivial conserved quantities

Numerical Implementation

• Tools: Julia package PauliStrings.jl (efficient computation of nested commutators using Pauli string representation)
• Translational Invariance: Store only unit cell rather than extensive operators
• Pauli Truncation Enhancement: Retain M Pauli strings with largest weights at each Lanczos step (Algorithm 1)

Comparison Methods

1. Matrix Product Operator (MPO) Simulation: As benchmark
2. Standard Dirichlet Boundary: Demonstrates non-physical reflection
3. Teretenkov et al. Method: Linear extrapolation + artificial dissipation (supplementary material comparison)

Evaluation Metrics

• Dynamical Symmetry Frequency Accuracy: Comparison with known results (XXZ) or analytical predictions (chaotic)
• Time Evolution Fidelity: TrO_0(0)O_0(t) vs exact/MPO results
• Locality Measure: Distribution of eigenvector in Krylov basis |ψ_0|² (darker indicates more local)

Experimental Results

Main Results

1. Ideal Case Verification (Figures 2-3)

Linear Growth b_n = n:

• Dynamics under open boundary conditions indistinguishable from infinite chain
• All eigenvalues satisfy Im(ω) = -2i (except two trivial roots ω = -i, -3i)
• Dirichlet boundary causes significant non-physical reflection at t≈5

Square-Root Growth b_n = √n:

• Open boundary conditions perfectly reproduce infinite chain dynamics
• Real part of spectrum nearly unchanged, imaginary part introduces small decay

2. XXZ Chain Results (Figures 4a,c)

• Dynamical Symmetry Recovery: Successfully extract known dynamical symmetry near ω = 12
  • L=26: Re(ω) ≈ 12.0, Im(ω) ≈ -0.5
  • Pauli truncation enhancement (L=50): convergence significantly improved
• Time Evolution: TrO_0(0)O_0(t) shows sustained oscillations (matches MPO to t≈3)
• Locality: Only real-valued dynamical symmetries are local (|ψ_0|² concentrated at small n)

3. Chaotic Chain Results (Figures 4b,d)

• No Eternal Symmetries: All modes decay, no purely real frequencies
• Decay Rate Clustering: Eigenvalues concentrated near Im(ω) ≈ -0.72i
• Agreement with Theory: -0.72 ≈ -2λ, where λ is Lanczos coefficient growth rate
• Exponential Decay: TrO_0(0)O_0(t) decays rapidly to zero (matches MPO)
• Non-locality: All modes highly non-local (|ψ_0|² uniformly distributed)

4. Steady-State Quench (Figure 5)

Quench from |ψ_0⟩ = 2^{-N/2}(|0⟩+|1⟩)^⊗N:

• Q_1 = ∑_i (s^+_i + s^-_i): rapid decay
• Q_5 (3rd order 5-local dynamical symmetry expansion): long-time oscillations
• Results from different L (24-30) converge well
• Standard truncation (gray line) exhibits non-physical behavior

Ablation Studies

Impact of Pauli Truncation (Supplementary Figure 4)

• Without Pauli Truncation: L_max = 26 (exact computation limit)
• With Pauli Truncation (retaining 2^22 strings): L_max = 50
  • Dynamical symmetry frequency: Re(ω) improves from ≈12.1 to ≈12.05
  • Im(ω) improves from ≈-0.8 to ≈-0.2

Convergence Analysis (Supplementary Figures 2-3)

• XXZ Model: Re(ω) stabilizes after L=10, Im(ω) continues improving
• Chaotic Chain: Spectrum essentially converges after L=30, deep-colored (large L) eigenvalues densely clustered near Im(ω)=-2λ line

Comparison with Teretenkov Method (Supplementary Figures 6-7)

Chaotic Chain:

• Teretenkov method (γ=0.2-0.4) eigenvalues do not converge to analytical result
• This method accurately reproduces clustering at Im(ω)=-2λ

XXZ Chain:

• Both methods show similar accuracy for local dynamical symmetries (ω=12)
• Predictions for non-local modes differ significantly

Experimental Findings

1. Quantitative Relationship Between Operator Growth and Thermalization: In linear growth (chaotic) case, decay rate precisely equals -2α, directly verifying theoretical predictions
2. Locality Criterion: |ψ_0|² in Krylov representation automatically provides quantitative measure of dynamical symmetry locality
3. Method Applicability:
   • Chaotic systems (linear growth): excellent performance, smoothness assumption well-satisfied
   • Integrable systems (square-root growth): can extract dynamical symmetries, but slower convergence
4. Computational Bottleneck: Primary limitation is computing high-order Lanczos coefficients (operators become highly non-local), Pauli truncation provides effective mitigation

Related Work

Krylov Space Methods

• Operator Complexity: Parker et al. 105 for probing operator growth in quantum chaos
• Hydrodynamics: Refs 97, 117-119 studying transport properties
• Floquet Systems: Refs 120-124 applications to periodically driven systems

Pseudo-mode Expansion

• Mori Projection: Early method for computing time correlation functions 92
• Recent Developments: Refs 96, 97 suggest adding dissipation + extrapolation or studying Ruelle-Pollicott resonances 98-104
• Limitations: Restricted to near-equilibrium, pseudo-locality not examined

Dynamical Symmetry Theory

• Eternal Equilibrium: Theoretical framework by Buča 82, state form ρ(t) = Z^{-1}exp(∑_u μ_u e^{iλ_u t}A_u)
• Spectrum-Generating Algebras: Refs 33, 34, 90, 91 related but not necessarily pseudo-local
• Time Crystals and Scars: Applications of dynamical symmetries in Refs 12, 29-44

Advantages of This Work

1. Parameter-Free: No introduction of artificial dissipation rates or other parameters
2. No Extrapolation Required: Direct use of known Lanczos coefficients
3. Pseudo-locality: Explicitly identifies relevant pseudo-local symmetries
4. Far from Equilibrium: Applicable to strong quenches and other non-equilibrium situations

Conclusions and Discussion

Main Conclusions

1. Method Effectiveness: Open boundary conditions successfully enable numerical computation of Liouvillian spectrum for infinite closed systems
2. Theoretical Insights:
   • Chaotic systems: linear Lanczos growth → all observables decay at rate -2α
   • Integrable systems: square-root growth → existence of eternal oscillating modes
3. Practical Value: Combined with Pauli string methods, provides powerful tool for computing many-body dynamics

Limitations

1. Smoothness Assumption: Method depends on smoothness of φ_n in Krylov representation
   • Chaotic systems: linear growth ensures smoothness
   • Other cases: verification required
2. Computational Complexity:
   • High-order operators O_n become extremely non-local
   • XXZ model can only be exactly computed to L=26
   • Pauli truncation introduces approximation but is effective
3. Long-Time Accuracy:
   • Incomplete convergence of imaginary part of dynamical symmetries causes long-time errors
   • Requires larger L or improved truncation strategies
4. Hot State Quench Difficulty: Quenching from thermal state ρ = e^{-βO_0} requires computing exponentials of non-commuting operator sums, numerically challenging (only steady-state handled)

Future Directions

1. Quantum Circuits: Generalization to cases where Liouvillian is upper Hessenberg matrix 107
2. Other Pseudo-local Symmetries:
   • Semi-local dynamical symmetries 137-139
   • Pseudo-local symmetries of quantum scars 30, 31, 36, 82
3. NMR Experimental Simulation: Direct computation of Fourier transform of Tr(Z_tot(t)Z_tot(0)) 140
4. Improved Convergence:
   • More refined truncation strategies
   • Adaptive Pauli truncation
   • Recursive algorithm optimization (supplementary material Algorithm)

In-Depth Evaluation

Strengths

1. Methodological Innovation ⭐⭐⭐⭐⭐

• Original Boundary Conditions: Based on natural physical assumption of wave function smoothness, distinct from artificial dissipation addition
• System-Environment Decomposition: Clever framework achieving local-nonlocal separation at operator space level
• Theory-Numerics Integration: Analytical solutions for ideal cases (Meixner/Hermite polynomials) complement numerical verification on actual systems

2. Theoretical Depth ⭐⭐⭐⭐⭐

• Operator Growth-Thermalization Connection: First quantitative proof of universal relation linear growth → -2α decay
• Pseudo-locality Quantification: Automatic identification of relevant dynamical symmetries through |ψ_0|²
• Eternal Equilibrium Framework: Deep integration with Buča's theory 82

3. Experimental Sufficiency ⭐⭐⭐⭐

• Multi-model Verification: Two typical cases—integrable (XXZ) and chaotic chains
• Multiple Comparisons: MPO, Dirichlet boundary, Teretenkov method
• Convergence Analysis: Detailed L-dependence and Pauli truncation studies
• Limitation: Lacks testing on more realistic physical systems (e.g., Hubbard model)

4. Technical Implementation ⭐⭐⭐⭐⭐

• Open-Source Code: Julia implementation provided (GitHub + PauliStrings.jl)
• Algorithm Optimization: Translational invariance, Pauli truncation, recursive algorithms and other practical techniques
• Reproducibility: Detailed supplementary material and parameter specifications

Weaknesses

1. Rigor of Smoothness Assumption

• Lack of Theoretical Proof: When is φ_n sufficiently smooth for linear approximation?
• Failure Cases: May fail under strong fragmentation or special boundary conditions
• Recommendation: Requires more rigorous mathematical proof or applicability criteria

2. Computational Scalability

• L Upper Limit: XXZ only reaches 26 (50 requires truncation), limiting high-precision applications
• System Size Dependence: N (system size) dependence not discussed
• 2D Generalization: Only 1D chains; higher dimensions challenging

3. Limited Experimental Verification

• Quench Types: Only steady-state quench, thermal state quench not implemented
• Time Scales: Long-time accuracy limited (t<5 for XXZ)
• Physical Quantities: Primarily operator auto-correlations; lacks other observables

4. Insufficient Comparison with Teretenkov Method

• Only in Supplementary Material: Main text lacks detailed discussion
• Parameter Sensitivity: Systematic study of γ and extrapolation schemes missing
• Advantage Analysis: When is this method superior? When equivalent?

Impact

Contribution to Field ⭐⭐⭐⭐⭐

• Paradigm Shift: Introduces open system concepts to closed system dynamics research
• Universal Tool: Applicable to broad range of quantum many-body problems
• Theoretical Foundation: Provides unified framework for understanding thermalization, time crystals, scars, etc.

Practical Value ⭐⭐⭐⭐

• Numerical Method: Directly applicable to computing actual system dynamics
• Experimental Guidance: Theoretical support for NMR and other experiments
• Limitation: Computational cost still restricts ultra-large system applications

Reproducibility ⭐⭐⭐⭐⭐

• Open-Source Code: Complete GitHub implementation
• Detailed Documentation: All algorithm details in supplementary material
• Parameter Transparency: All experimental settings clearly specified

Applicable Scenarios

Ideal Scenarios

1. 1D Local Interaction Systems: Spin chains, lattice bosons, etc.
2. Chaotic Systems: Linear Lanczos growth ensures method effectiveness
3. Short to Intermediate Dynamics: t<10τ (characteristic time scale)
4. Local Observables: Measurement of k-local operators

Scenarios Requiring Caution

1. Strongly Fragmented Systems: Smoothness assumption may fail
2. Long-Range Interactions: Krylov space structure may differ
3. Extremely Long-Time Dynamics: Incomplete imaginary part convergence causes error accumulation
4. High-Dimensional Systems: Computational cost increases dramatically

Inapplicable Scenarios

1. Open Systems: Lindbladian framework more appropriate
2. Strong Measurement Systems: Requires different theoretical tools
3. Classical Limit: Overcomplicated; classical methods simpler

Technical Details Supplement

Key Formula Derivations

Master Equation Form of Boundary Conditions (Supplementary Material S20)

Exact time evolution Ȯ = iH,O, hard boundary condition corresponds to:

Ȯ = iH,O - 1/(2N) Tr(O_L† iH,O)O_L

Substituting open boundary condition:

Ȯ = iH,O - 1/(2N) Tr(O_L† iH,O)O_L + ib_{L+1}Tr(O_† O)O_L - 2b_{L+1}Tr(O_L† O)O_L

Demonstrates how boundary condition balances operator flow into environment

Meixner Polynomial Solution (Supplementary Material S3-S7)

For linear growth b_n = n, eigenvector:

φ_n(ω) = i^n M_n(ω)/n!

where M_n satisfies recursion: M_{n+1}(x) = xM_n(x) - n²M_(x)

Boundary condition yields frequency constraint: (ω + 2i(L+1))M_L(ω) = (2L+1)LM_(ω)

Trivial roots ω = -i, -3i correspond to φ_n ∝ 1 and φ_n ∝ 2n-1

Numerical Techniques

Recursive Algorithm (Supplementary Material Algorithm)

Improved dynamical symmetry extraction:

1. Randomly initialize k-local O_0
2. Diagonalize L̃, select eigenmode with frequency ≈ω as A_ω
3. Truncate A_ω to k-local, update O_0 ← truncate(A_ω)
4. Repeat until convergence

XXZ model with k=3 converges in 4 steps: Re(ω) from 12.14→12.02

Selected References

Methodological Foundations

• 82 B. Buča, PRX 13, 031013 (2023) - Eternal equilibrium theory
• 105 D. Parker et al., PRX 9, 041017 (2019) - Krylov complexity
• 96,97 A. Teretenkov, O. Lychkovskiy et al., PRB (2024) - Pseudo-mode expansion

Application Domains

• 12 B. Buča et al., Nat. Commun. 10, 1730 (2019) - Time crystals
• 29-44 Quantum scars literature series
• 119 N. Loizeau et al., SciPost (2025) - PauliStrings.jl

Theoretical Background

• 1 L. D'Alessio et al., Adv. Phys. 65, 239 (2016) - Thermalization review
• 86-89 Pseudo-locality theory (Prosen, Doyon, et al.)

Overall Assessment: This is an excellent paper with important innovations in quantum many-body dynamics (rating 4.5/5). The core method (open boundary Krylov chain) is elegant and physically intuitive, theoretical derivations rigorous (analytical solutions for ideal cases), and numerical verification comprehensive. Main contributions establish quantitative connection between operator growth and thermalization, providing practical numerical tools. Primary limitations concern applicability range of smoothness assumption and computational scalability. For researchers studying non-equilibrium dynamics, time crystals, quantum scars, and related frontier topics in quantum many-body systems, this is essential reading.