Scalable simulation of non-Markovian quantum transport by stochastic-phase bath reduction

AI summary

Overview

Research area: Quantum physics — non-Markovian open quantum dynamics and computational simulation of quasiparticle transport in molecular materials.

Technical level: Advanced. The paper combines a rigorous operator-theoretic error proof with large-scale numerical benchmarks and assumes familiarity with Holstein models, bath correlation functions, trace norm, and Lindblad dynamics.

Scope (one sentence): The paper introduces the Stochastic Phase Algorithm (SPA), which replaces many independent local vibronic baths with a few shared baths carrying random site-dependent phases, and demonstrates system-size-independent error scaling plus benchmarks on lattices up to 100 × 100 and molecular systems including BChl aggregates and rubrene.

What This Paper Is About

Simulating non-Markovian quantum transport in extended molecular systems is difficult because every electronically active site is coupled to its own structured vibronic environment, so the Hilbert-space dimension grows exponentially with system size. Existing systematically improvable methods (influence functionals, hierarchies, tensor networks) become expensive as environmental memory and system–bath entanglement grow, while semiclassical methods scale better but do not retain system–bath entanglement or converge systematically to exact quantum dynamics. The paper's goal is to reduce the spatial number of baths without discarding quantum system–bath dynamics, and to prove that the resulting error is controlled independently of system size.

Key Contributions

1. The Stochastic Phase Algorithm (SPA). It replaces the $N_{\text{s}}$ independent local Gaussian baths of a single-quasiparticle Holstein model with $R$ shared quantum baths dressed by static, stochastic site-dependent phases, each realization retaining the local bath spectrum and supporting system–bath entanglement.

2. An exactness statement for two-point correlations plus a size-independent error bound. The authors prove that phase averaging reproduces the full two-point bath-correlation matrix for every $R$, and that higher-order cross-site correlation errors are suppressed as $1/R$. They derive the fixed-time trace-norm bound $|\hat{\rho}{\text{s}}(t)-\hat{\rho}{\text{s}}^{(R)}(t)|{\text{tr}} \le A_t/R$, where $A_t$ depends on time and the bath correlation functions but not on $N{\text{s}}$ or network connectivity; for a bounded observable the error is at most $|\hat{O}|_\infty A_t/R$.

3. Combination with coupled-Lindblad spectral compression. SPA reduces the number of spatial bath copies, while coupled-Lindblad compression reduces the number of modes inside each shared bath by fitting damped quantum modes to the bath correlation function over the simulation window, giving state dimensions of $N_{\text{s}}d^{nR}$ (unitary shared baths) and $N_{\text{s}}d^{n_{\text{b}}R}$ (after spectral compression) instead of $N_{\text{s}}d^{nN_{\text{s}}}$.

4. Benchmarks and open software. SPA is validated against exact unitary dynamics on a 3 × 3 lattice, run as a full state-vector quantum-bath simulation on a 100 × 100 lattice, and compared with modular path-integral and TD-DMRG benchmarks for BChl aggregates and rubrene. The authors release the GPU-parallel scalabath Python package and use the RealTimeBath package for coupled-Lindblad fits.

Main Findings

• Error scaling: Phase averaging restores the two-point bath correlation functions exactly, since $\mathbb{E}\theta[R^{-1}\sum_a r_i^{(a)}r_j^{(a)*}]=\delta{ij}$. Fourth-order cross-site contributions are controlled by $\mathbb{E}\theta[G{ij}^{(R)}G_{ji}^{(R)}]=1/R$, and the overall trace-norm error is bounded by $A_t/R$ with a prefactor independent of system size and connectivity. At fixed time and target accuracy, enlarging the network does not require more bath copies.

• Two-point correlations alone are insufficient: the phase-averaged Gaussian baths need not remain Gaussian, so the proof accounts for all higher-order errors rather than only matching the two-point correlation matrix.

• 3 × 3 lattice convergence: With nearest-neighbor hopping $V=100\ \mathrm{cm}^{-1}$, mode frequency $\omega_1=1000\ \mathrm{cm}^{-1}$, $T=300\ \mathrm{K}$, and an initial state localized at the central site, the SPA error in central-site population and mean-squared displacement grows with $\lambda$ but decreases systematically with $R$. TCL2 is comparably accurate only at weaker coupling; for $\lambda/V>2$ it breaks down and produces nonphysical dynamics that violate positivity, whereas the completely positive SPA remains robust through the intermediate-coupling regime.

• Scalability to 100 × 100: With disordered site energies sampled uniformly in $[-W,W]$ with $W=200\ \mathrm{cm}^{-1}$ and $\lambda/V=8$, SPA with $R=1$ qualitatively agrees with a renormalized tight-binding reference using the Lang–Firsov hopping $V'=V\exp[-g_1^2\coth(\beta\omega_1/2)]$; both show slow diffusion over 10 ps due to strong-coupling suppression and disorder frustration. The paper notes that state-of-the-art tensor-network simulations retaining site-resolved environments have demonstrated long-time dynamics only for tens of sites.

• Exciton migration in a photosynthetic complex: For homogeneous and staggered 19-site B850-chain models with 50 vibrational modes per site, $\lambda=217.66\ \mathrm{cm}^{-1}$, and $T=300\ \mathrm{K}$, the 50-mode local bath is compressed into a shared six-mode dissipative bath with $\tau=100\ \mathrm{fs}$. SPA with $R=1$ closely reproduces modular path-integral results, capturing coherent spreading and vibrational damping in the homogeneous chain and energy-biased redistribution in the staggered chain; the bath-compression error is negligible over this window, and small deviations beyond $t\approx30\ \mathrm{fs}$ are consistent with the finite-$R$ error.

• Charge transport in rubrene: For the homogeneous 1D Holstein model with $V=83\ \mathrm{meV}$, nine vibrational modes per site, $\lambda=73\ \mathrm{meV}$, and $\lambda/V\simeq0.88$, SPA directly simulates a chain of $L=200$ sites and computes the $b$-axis mobility. It closely reproduces the TD-DMRG mobility from 200 to 400 K, whereas Ehrenfest dynamics shows systematic deviations and data-informed quantum-classical dynamics shows nonmonotonic residual errors. The authors note that their Ehrenfest and DIQCD setups differ from those in the earlier reference, which used shorter propagation times, smaller bosonic cutoffs, and $D=7\ \mathrm{\AA}$.

• Physical consistency: Phase averaging preserves complete positivity, and the unitary thermal-bath construction recovers detailed balance in the Davies limit via the thermal KMS relation.

• Stated limitation: The size-independent bound generally does not extend to many-body settings where the relevant electronic sector grows with system size.

Methodology in Plain English

The starting model is a tight-binding network of $N_{\text{s}}$ sites with arbitrary connectivity, where each site carries its own set of $n$ bosonic modes with local diagonal coupling. Instead of keeping all $N_{\text{s}}$ separate environments, SPA keeps only $R$ copies of the bath and couples each copy to every site, multiplying each site–bath coupling by a random phase factor $r_i^{(a)}=e^{i\theta_i^{(a)}}$ and scaling the coupling by $1/\sqrt{R}$. The phases are drawn independently and uniformly from $[0,2\pi)$.

For each fixed set of phases, the dynamics is ordinary unitary quantum evolution of the system plus the shared baths, started from a thermal bath state. One then averages the reduced system density operator over many phase realizations. Because the phase factors cancel on average when site indices match and vanish when they differ, the averaged model reproduces the diagonal bath correlation matrix of the original independent baths. Randomness at fourth order and beyond does not cancel exactly, but it shrinks as more shared baths are used, which is what the $1/R$ bound quantifies.

For baths with many modes or continuous spectral densities, the authors further replace each shared bath by a few damped quantum modes that exchange excitations, fitting a Hermitian mode-coupling matrix, a positive damping matrix, and a coupling vector so that the model correlation function matches the target over the simulation window. The auxiliary modes start in vacuum, with temperature encoded in those fitted matrices. The resulting dynamics is a Lindblad equation, simulated with a Monte Carlo wavefunction propagator averaged jointly over quantum trajectories and static phases.

The numerical work compares SPA against exact unitary propagation where feasible (3 × 3), against a renormalized tight-binding reference at strong coupling (100 × 100), against modular path-integral results (19-site BChl chains), and against TD-DMRG, Fermi's golden rule, Boltzmann theory, Ehrenfest dynamics, and DIQCD (rubrene).

Why This Matters

The work removes a structural bottleneck in non-Markovian transport simulation: the need for an independent quantum bath per site. Because the error bound's prefactor does not scale with system size or connectivity, larger networks can be simulated at fixed accuracy and fixed computational cost in the bath sector. This lets quantum-bath dynamics reach mesoscale problems previously handled mainly by semiclassical methods or by quantum-bath approaches whose favorable scaling relies on quasiparticle localization.

Applications discussed in the paper (some as demonstrated benchmarks, some as proposed extensions):

• Organic electronics: Predicting charge-carrier mobility in high-mobility organic semiconductors, demonstrated for the $b$-axis of rubrene in the intermediate-coupling regime where weak- and strong-coupling theories fail.
• Photosynthetic light harvesting: Modeling exciton migration in antenna complexes such as LH2, with 9 BChl $a$ molecules in the B800 ring and 18 in the B850 ring.
• Molecular-qubit arrays: Proposed extension to spin-phonon coupling in molecular qubits.
• Broader open-system modeling: Proposed extension to fermionic environments and nonlocal electron-vibronic coupling, including Holstein–Peierls models, through correlated stochastic phases and auxiliary-bath couplings.

Industry relevance: The released GPU-parallel scalabath package, which supports systems with thousands of molecular sites, targets computational materials modeling where predictive transport properties of molecular semiconductors and molecular-scale quantum hardware are needed at scales beyond exact methods.

Future Directions

• Multidimensional transport: The authors state that organic systems can exhibit strongly anisotropic tunneling requiring 2D or 3D transport models, and that such multidimensional transport will be explored in future work.
• Tensor-network representations: Adopting a tensor-network representation of the system–bath state is suggested as a route to further scalability beyond the full state-vector representation.
• Broader environment types: Extending the framework to fermionic environments and to nonlocal electron-vibronic coupling, including Holstein–Peierls models, via appropriately correlated stochastic phases and auxiliary-bath couplings.
• Many-body regimes: The size-independent error bound is explicitly noted not to extend generally to many-body settings where the electronic sector grows with system size, leaving that regime open.
• Mode-count scaling: Standard unitary discretizations of continuous baths grow linearly in mode count with simulation duration, whereas coupled-Lindblad compression can achieve polylogarithmic growth; exploiting this fully for long-time mesoscale simulations is a natural next step.

Target Audience

Advanced researchers and graduate students in open quantum systems, computational physical chemistry and condensed-matter physics, and organic-electronics materials modeling. It is most valuable to method developers working on non-Markovian quantum dynamics, users of path-integral, hierarchy, TD-DMRG, or tensor-network approaches who need larger system sizes, and computational scientists interested in practical GPU implementations of quantum-bath transport simulation. Readers without a background in bath correlation functions, trace-norm error analysis, and Holstein-type models will find the theoretical sections demanding.

Authors’ abstract

Accurate simulations of non-Markovian quantum transport remain limited to small systems, and a key bottleneck is that each system site requires an independent set of quantum bath orbitals. We introduce the Stochastic Phase Algorithm (SPA), which replaces the independent local Gaussian baths in single-quasiparticle Holstein models by $R$ shared quantum baths with stochastic site-dependent phases. We prove that phase averaging reproduces the target two-point bath correlation matrix for every $R$, while increasing $R$ systematically suppresses errors from higher-order cross-site correlations. We obtain a fixed-time $O(R^{-1})$ trace-norm error bound with a prefactor independent of system size and connectivity. This spatial compression applies directly to unitary baths and can be combined with coupled-Lindblad spectral compression [Phys. Rev. Lett. 136, 090403 (2026)], which represents each shared environment by a very small number of coupled, damped quantum modes fitted to its thermal bath correlation over the simulation window. Numerically, SPA converges against direct unitary benchmarks on a $3\times3$ lattice, enables full-state-vector quantum-bath dynamics on a $100\times100$ lattice, and, with $R=1$, closely reproduces near exact benchmarks for exciton populations in bacteriochlorophyll aggregates and carrier mobility in rubrene.