Numerically Exact Floquet Propagator

• The paper presents a method using high-order Magnus and Chebyshev expansions to construct Floquet propagators with machine-precision control of errors.
• It details algorithms for closed systems via time-slicing and parallelization, and for open systems using MPO-based Trotterization to capture non-Markovian dynamics.
• The approach enables precise extraction of quasienergies, Floquet states, and time-resolved observables, underpinning analysis of energy flows and entanglement in quantum systems.

A numerically exact Floquet propagator is a computational construct that enables the precise simulation of periodically driven quantum systems, both isolated and dissipative, over discrete stroboscopic times. In the context of closed systems, the Floquet propagator is a unitary operator—constructed by integrating the Schrödinger equation over one driving period—that encapsulates the full periodic evolution, allowing for the extraction of quasienergies and Floquet states. For non-Markovian open quantum systems, an analogous construct arises at the level of the reduced density matrix, formulated as a matrix product operator (MPO) acting as a discrete-time superoperator, and constructed directly from the path-integral influence functional. These numerically exact schemes are designed to control all sources of error—time-discretization, operator expansion, and numerical round-off—to the level of machine precision, and are applicable in regimes where analytic or perturbative techniques fail.

1. Theoretical Foundation of Floquet Propagators

A Floquet propagator describes the time evolution over one period $T$ of a quantum system governed by a $T$-periodic Hamiltonian $H(t)$. For closed systems, this propagator is given by

$$U_F = \mathcal{T} \exp\left[ -\frac{i}{\hbar} \int_0^T H(\tau) d\tau \right],$$

where $\mathcal{T}$ denotes time-ordering. The eigenvalues $e^{-i\theta_\alpha}$ and eigenvectors $|\varphi_\alpha\rangle$ of $U_F$ yield the quasienergies $\epsilon_\alpha = \hbar \theta_\alpha/T$ and the Floquet states, which encode the stroboscopic dynamics.

For open quantum systems with periodic driving and strong non-Markovian dissipation, the system’s reduced evolution over one period is described not by a unitary, but by a superoperator Floquet map acting on the system density matrix. This map is derived from a path-integral-based influence functional, which is compressed into a periodic MPO and then contracted over one period to yield the Floquet propagator $\mathcal{U}_F$ (Mickiewicz et al., 11 Nov 2025).

2. Numerically Exact Construction for Closed Systems

The numerically exact procedure for constructing $U_F$ in finite-dimensional, closed quantum systems proceeds as follows (Laptyeva et al., 2015):

1. Time-Slicing and Magnus Propagator: Partition the period $T$ into $M$ slices of length $h=T/M$, defining time nodes $t_j = j h$. For each slice $[t_{j-1}, t_j]$, approximate the exact propagator via a single exponential using a truncated Magnus expansion:

$$U_j \approx \exp\left[-\frac{i}{\hbar} \Omega_j\right],$$

where $\Omega_j$ is calculated up to sixth order for optimal balance between cost and accuracy, e.g., via the Blanes–Casas commutator-free formula which achieves $\mathcal{O}(h^7)$ accuracy.

1. Chebyshev Expansion of Exponentials: The matrix exponential $\exp(-i\Omega_j/\hbar)$ is computed efficiently via a Chebyshev polynomial expansion after spectral rescaling (so that $\text{spec}(\tilde{\Omega}_j) \subset [-1,1]$). Expansion to order $L \approx R + 10$–$20$, with $R = (\Delta E \cdot h)/\hbar$, ensures convergence to machine precision $\varepsilon \sim 10^{-14}$.
2. Assembly of the Floquet Operator: The full Floquet operator is constructed as $U_F = U_M U_{M-1} \cdots U_1$ by repeated propagation of the identity matrix, leveraging block-wise parallelization across $P$ nodes (MPI), with BLAS-accelerated matrix-matrix multiplications.
3. Diagonalization and Floquet State Reconstruction: The completed $U_F$ is diagonalized (single-call cost $\mathcal{O}(N^3)$), e.g., with threaded LAPACK/MKL routines for $N$ up to $10^4$. The complete time-dependent Floquet states $|\Phi_\alpha(t)\rangle$ are reconstructed by re-propagating the eigenvectors through all $M$ time steps, enabling detailed observable and spectral analyses.

This workflow achieves strictly unitary evolution (modulo numerical round-off), exploits highly parallelizable step-wise structure, and yields phase-overlap errors $\Sigma_\alpha \lesssim 10^{-12}$ for $M \sim 10^4$ (i.e., $\Sigma_\alpha \propto M^{-6}$ for 6th order Magnus–Chebyshev expansion) (Laptyeva et al., 2015).

3. Numerically Exact Construction for Open Systems: Periodic Influence MPO

For open systems described by a system-plus-bath model (e.g., spin-boson problems), numerically exact Floquet propagation is achieved as follows (Mickiewicz et al., 11 Nov 2025):

1. Feynman–Vernon Functional Discretization: The full reduced dynamics is formulated via the path integral, involving an influence functional $\mathcal{F}[s^+, s^-]$ encoding the effect of the environment. Discretizing the time axis and exploiting underlying periodicity, the influence tensor is cast as a periodic MPO,

$$\mathcal{I}_{s^+_N s^-_N, \ldots, s^+_1 s^-_1} = \sum_{\{r_k\}} W^{[N]\,s^+_N, s^-_N}_{r_{N-1}, r_N} \cdots W^{[1]\,s^+_1, s^-_1}_{r_0, r_1},$$

where $W^{[n+M]} = W^{[n]}$. Bond indices $r_n$ scale as MPO bond-dimension $\chi \sim 200$–$400$ for practical bath spectra.

1. Embedding the Driven System via Trotterization: Incorporation of the driven system dynamics is effected via symmetric Trotter decomposition. Each block $Q_n^{\mu\mu'}$ is constructed by interleaving exact local system propagators (solved numerically) with the bath-mediated uniform $q^{(\mu,\mu')}$ blocks, so as to assemble the time-dependent $Q_n$. Periodicity ensures $Q_{n+M} = Q_n$.
2. Formation and Diagonalization of Floquet Superoperator: The one-period Floquet map is

$Q_F = Q_M Q_{M-1} \cdots Q_1,$

acting as a superoperator on the system density matrix. Its eigenspectrum is obtained via iterative MPO eigensolvers (Arnoldi/Lanczos), allowing access to the asymptotic periodic state ($|\lambda_1| = 1$), decay rates ($|\lambda_{2,3,\dots}| < 1$), and stroboscopic evolution,

$\rho(kT) = \sum_i c_i \lambda_i^k |\phi_i\rrangle.$

1. Numerical Control and Cost: The time-step $\delta t$ is chosen so that Trotter errors $\mathcal{O}((\delta t)^3)$ are negligible, and the MPO bond dimension $\chi$ systematically increased until two-point correlators and local populations converge within a prescribed threshold ($10^{-6}$). Contraction and diagonalization, scaling as $\mathcal{O}(M d^2 \chi^3)$ and $\mathcal{O}(k_{\text{it}} d^2 \chi^3)$ respectively, are computationally efficient for $d$ up to several tens and $\chi \lesssim 400$ (Mickiewicz et al., 11 Nov 2025).

4. Parallelization Strategy and Computational Scaling

The independence of the sub-propagators in both approaches makes the numerically exact Floquet propagator highly amenable to parallelization:

• Closed Systems:

The identity operator is partitioned into $P$ blocks, each propagated in parallel. MPI distributes block workloads across nodes, while within each node OpenMP/MKL maximizes efficiency of BLAS operations. Memory per node is $2$–$3 N^2/P$ complex doubles, sufficient for storing block and propagator data. Reported timings for $M = 10^2, N_{\text{Ch}} = L = 50$ demonstrate sub-hour runs for $N \le 10^3$, with multi-hour scalability to $N = 10^4$. Strong scaling is observed, with efficiency increasing with larger problem size ($N$) and time-slice number ($M$), owing to improved BLAS utilization (Laptyeva et al., 2015).

• Open Systems (MPO-Based):

Computational cost for an individual $Q_n$ block scales as $\mathcal{O}(d^2 \chi^2)$, with full contraction and eigendecomposition costs as $\mathcal{O}(M d^2 \chi^3)$ and $\mathcal{O}(k_{\text{it}} d^2 \chi^3)$ respectively. The fixed MPO structure allows for efficient contraction on distributed memory platforms, and iteration counts for eigensolvers are modest ($k_{\text{it}} \sim 20$–$50$).

A practical implication is that exact Floquet analysis of systems with $N \sim 10^4$ (closed) or Liouville space dimension $d$ up to $10^2$ (open) is computationally accessible on modest clusters, with parallel scaling improving as system size and periodic slicing increase.

5. Applications to Quantum Dynamics and Observables

The numerically exact Floquet propagator framework enables precise stroboscopic and intra-period evaluation of:

• Quasienergies and Floquet States:

Direct computation of the Floquet spectrum and eigenmodes, quantifying resonance structure, avoided crossings, and topological signatures in time-periodic systems (Laptyeva et al., 2015).

• Time-Resolved Observables:

Arbitrary operator expectation values can be evaluated as

$$\langle A \rangle_\alpha = \frac{1}{T}\int_0^T \langle \Phi_\alpha(t) | A(t) | \Phi_\alpha(t) \rangle dt$$

in closed systems, or as multi-time correlation functions in open systems, via operator insertions at specific MPO sites.

• Heat Currents in Driven Spin-Boson Models:

In the driven spin–boson model, the Floquet-MPO method allows accurate computation of period-averaged heat currents,

$$\overline{j}(\omega) = 2J(\omega)\omega \int_0^\infty d\tau\, \mathrm{Im} \left\{ [1+2n_B(\omega)]\sin(\omega\tau) + i\cos(\omega\tau) \right\} \overline{C}(\tau),$$

with $\overline{C}(\tau)$ the period-averaged correlator. Numerical findings reveal qualitative distinctions between longitudinal and transversal driving, such as broad heat absorption bands and sharply peaked spectral features, with total dissipated power $\overline{I}$ exhibiting explicit non-equilibrium behavior (Mickiewicz et al., 11 Nov 2025).

• Floquet-Engineered Entanglement:

For two qubits in a common bath, periodic driving stabilizes a non-trivial steady-state concurrence, with maximal entanglement arising at resonance $\omega_d \approx 2 \Omega$, $\epsilon_d \approx 1.1 \Omega$. This phenomenon is directly attributable to a long-lived excited eigenmode of the Floquet superoperator, showcasing the control of steady states beyond reach for purely Markovian or perturbative approaches.

6. Accuracy, Limitations, and Error Analysis

Numerical exactness is ensured by:

• Systematic Control of Truncation Errors:

For closed systems, increasing $M$ and the order of the Magnus and Chebyshev expansions reduces errors below $10^{-12}$, with strict unitarity maintained up to round-off.

• Convergence in MPO Bond Dimension and Timestep:

In open systems, $\chi$ and $\delta t$ are systematically increased/decreased until physical observables stabilize within target tolerances, e.g., $10^{-6}$ for local observables, allowing full control over non-Markovian dynamics.

No evidence is provided for fundamental scaling barriers up to $N\sim 10^4$ (closed), or $d\sim 100$ and $\chi\sim 400$ (open), barring hardware limitations. Iterative eigensolvers and block-parallelization mitigate the cubic scaling bottlenecks traditionally facing large-scale Floquet analysis (Laptyeva et al., 2015, Mickiewicz et al., 11 Nov 2025).

7. Significance and Broader Context

The construction of numerically exact Floquet propagators, whether as unitary matrices (closed systems) or as periodic MPOs (strongly damped open systems), represents a definitive tool for the high-fidelity simulation of quantum dynamics under periodic driving. These frameworks extend rigorous time-dependent simulation capability to regimes of strong driving, strong dissipation, and nontrivial system-bath coupling that are inaccessible to standard Markovian or perturbative master equations.

Applications range from spectral and entanglement engineering to the computation of heat currents and steady-state properties, with demonstrated accuracy and scalability. Consequently, these approaches enable the systematic investigation of quasienergy spectra, dynamical localization, emergent entanglement mechanisms, and non-equilibrium energy flows in both fundamental and applied contexts across condensed matter, quantum optics, and quantum information theory.

The systematic control over all approximation errors inherent in these propagators underpins their "numerically exact" designation, underscoring their relevance for the computational paper of driven quantum matter and dissipative quantum systems at arbitrary coupling and periodicity scales (Laptyeva et al., 2015, Mickiewicz et al., 11 Nov 2025).