- The paper presents a formal treatment showing how a closed SYK system thermalizes following a non-equilibrium quench, revealing parameter-dependent thermal states.
- It employs Schwinger-Dyson and Kadanoff-Baym equations coupled with a causal-stepping predictor-corrector integrator to solve complex integro-differential dynamics.
- Key outcomes include the extraction of final thermal temperatures and exponential relaxation rates, underscoring non-universal thermalization behavior in quantum chaotic systems.
Thermalization Dynamics of the Closed Majorana SYK Model in the Thermodynamic Limit
The paper addresses the real-time quantum relaxation dynamics of a large-q closed Majorana Sachdev-Ye-Kitaev (SYK) system under a non-equilibrium quench. The SYK model, with its lack of quasiparticle description, maximal chaos, and closed-form solvability in the large-N limit, serves as a valuable setting to probe the fundamentals of quantum thermalization in isolated, non-integrable many-body systems. The model under study incorporates two SYK-type terms: a q-body all-to-all interaction (the "interaction term") and a random quadratic hopping ("kinetic" term). The system is initialized in equilibrium with the interaction term, and subsequently, at t=0, a kinetic term of tunable strength is instantaneously switched on, resulting in a mixed post-quench Hamiltonian.
In this architecture, the two terms compete, and the model remains exactly treatable in the N→∞, q→∞ limits via closed Schwinger-Dyson and Kadanoff-Baym equations. The setup allows for parameterization of the quench via the initial temperature and the kinetic term strength J2, offering access to a broad regime from weak to strong quantum quenches.
Figure 1: Sketch of the Schwinger-Keldysh contour C=C++C−+Cimag employed in the real-time non-equilibrium path integral formalism.
Analytical Framework: Schwinger-Dyson and Kadanoff-Baym Equations
The central analytical development uses the Keldysh non-equilibrium path-integral formalism. Starting from the disorder-averaged partition function, the authors derive closed dynamical equations for the Green's functions, exactly valid in the thermodynamic (N→∞) and large-q limits. The resultant Kadanoff-Baym equations are nonlinear, nonlocal, and non-Markovian, relating greater and lesser Majorana Green’s functions, and include explicit dependence on the full quantum history along the Keldysh contour.
To treat the large q regime analytically and numerically, an exponential ansatz for the Green’s functions is employed, which collapses the computational complexity and yields a tractable Volterra-type partial integro-differential equation for a complex scalar function g(t1,t2). This formulation admits analytic treatment in two limiting cases: (1) pure SYK interaction—leading to analytic thermal Green’s functions with instantaneous relaxation, and (2) quench to a pure kinetic (hopping) Hamiltonian—yielding instantaneous thermalization to an infinite-temperature state, irrespective of the initial conditions.
Figure 2: Sketch of the fan-like causal propagation and stepping scheme in the two-time (t1,t2) plane employed to numerically integrate the Kadanoff-Baym equations.
Numerical Solution: Causal Stepping and Predictor-Corrector Integrator
The general mixed quench scenario, involving both interaction and kinetic terms post-quench, defies full analytic treatment due to the complex memory structure and strongly non-Markovian character of the Volterra equations. The authors develop a dedicated numerical solver based on a causal-stepping predictor-corrector integration scheme in the two-time plane. By exploiting the underlying symmetry and structure of the equations (Majorana symmetry, conjugation, causal propagation), the computational cost is reduced and precise energy conservation is achieved to sub-percent accuracy throughout the time evolution.
The predictive-corrective approach is benchmarked against the analytically solvable limiting cases—both in terms of the time-dependent energy and norm-wise Green’s function errors—and achieves excellent agreement within a few tenths of a percent across the temperature, quench, and coupling phase space.



Figure 3: Predictor-corrector performance in the limiting cases: numerical versus analytic energies and Green’s function errors for equilibrium and kinetic quench protocols.
Non-Equilibrium Relaxation and Emergence of Thermal Stationarity
Time-evolution visualizations and analysis reveal clear evidence of thermalization at late times for the general mixed quench. After initial transient non-equilibrium behavior following the quench, the two-time Green’s function loses dependence on the average time T=(t1+t2)/2 and attains stationarity in the relative time t=t1−t2; i.e., g(T,t)→gβf(t) for T→∞. This is verified quantitatively by demonstrating vanishing derivatives with respect to T at fixed t and that the stationary solution satisfies the equilibrium Keldysh equations at a well-defined effective inverse temperature βf.
The stationary Green’s function is used to extract observables, most notably the final thermal energy, and compared directly to the expectation value computed from an imaginary-time equilibrium solution at the numerically fitted βf. This approach yields an accurate determination of the effective temperature as a function of quench strength and initial state. The stationary solution is unique and independent of trajectory details for fixed parameters, in agreement with the Picard-Lindelöf uniqueness theorem for initial value ODEs.
Quantitative Results: Final Temperature and Thermalization Rate
The main numerical results characterize the dependence of both the effective post-quench thermal temperature and the thermalization rate on the initial temperature and the quench parameter J2. Strong numerical convergence is achieved:
Theoretical and Practical Implications
This work establishes definitively that a non-integrable, maximally chaotic, all-to-all interacting quantum system—that remains exactly solvable at N→∞, q→∞—exhibits thermalization under nontrivial quantum quenches, but the resulting post-quench state is strongly non-universal and only approaches instantaneous relaxation in specific limiting cases (pure interaction or pure hopping in the final Hamiltonian). The dynamics reveal that the route to thermalization is sensitive to both the degree of the quench and the microscopic model parameters, even in a system with maximal quantum chaos and no quasiparticle description.
Practical implications are substantial for the theoretical description of quantum ergodicity and for quantum simulation architectures that attempt to realize or harness such models. The strong dependence of thermalization time on perturbation and temperature implies that generic quantum simulators of non-integrable models will see energy absorption and thermalization timescales that are highly tunable and potentially observable. The numerical methods developed—particularly causal-stepping two-time integrators for non-Markovian IDEs—provide tools readily adaptable to a range of quantum many-body and open system settings.
Figure 5: Goodness of exponential fit for kinetic energy relaxation, validating the extraction of the thermalization rate γ across parameter space.
Conclusions
The work rigorously demonstrates finite-time, parameter-dependent thermalization in quantum-coherent, closed SYK systems at the analytic level, quantified via non-equilibrium real-time dynamics. The critical results are the explicit computation of final temperatures and relaxation rates as a function of quench protocol and initial conditions, verification of unique thermal stationary states, and numerical evidence of rich thermalization phenomenology. The methodology and approach are directly extensible to complex SYK networks, higher dimensions, and other non-equilibrium quantum dot and lattice models. This study thereby contributes a robust analytic and computational foundation for understanding chaos-driven thermalization in strongly correlated quantum matter.