Static disorder plays a crucial role in the electronic dynamics and spectroscopy of complex molecular systems. Traditionally, obtaining observables averaged over static disorder requires thousands of realizations via direct sampling of the disorder distribution, leading to high computational costs. In this work, we extend the auxiliary degree-of-freedom based matrix product state (MPS) method to handle system-bath correlated thermal equilibrium initial states, which can capture static disorder effects using a one-shot quantum dynamical simulation. We validate the effectiveness of the extended method by computing the dipole–dipole time correlation function of the Holstein model relevant to the emission spectrum of molecular aggregates. Our results show that the one-shot method is very accurate with only a moderate increase in MPS bond dimension, thereby significantly reducing computational cost. Moreover, it enables the generation of a much larger number of samples than the conventional direct sampling method at negligible additional cost, thus reducing statistical errors. This method provides a broadly useful tool for calculating equilibrium time correlation functions in system-bath coupled models with static disorder.
Static disorder plays a crucial role in the electronic dynamics and spectroscopy of complex molecular systems. Traditionally, obtaining observables averaged over static disorder requires thousands of realizations via direct sampling of the disorder distribution, leading to high computational costs. In this work, we extend the auxiliary degree-of-freedom based matrix product state (MPS) method to handle system-bath correlated thermal equilibrium initial states, which can capture static disorder effects using a one-shot quantum dynamical simulation. We validate the effectiveness of the extended method by computing the dipole–dipole time correlation function of the Holstein model relevant to the emission spectrum of molecular aggregates. Our results show that the one-shot method is very accurate with only a moderate increase in MPS bond dimension, thereby significantly reducing computational cost. Moreover, it enables the generation of a much larger number of samples than the conventional direct sampling method at negligible additional cost, thus reducing statistical errors. This method provides a broadly useful tool for calculating equilibrium time correlation functions in system-bath coupled models with static disorder.
Understanding the emergent system-bath correlations in non-Markovian and non-perturbative open systems is a theoretical challenge that has benefited greatly from the application of Matrix Product State (MPS) methods. Here, we propose an autonmously adapative variant of the one-site Time-Dependent-Variational-Principle (1TDVP) method for many-body MPS wave-functions in which the local bond-dimensions can evolve to capture growing entanglement 'on the fly'. We achieve this by efficiently examining the effect of increasing each MPS bond-dimension in advance of each dynamic timestep, resulting in an MPS that can dynamically and inhomogeneously restructure itself as the complexity of the dynamics grows across time and space. This naturally leads to more efficient simulations, oviates the need for multiple convergence runs, and, as we demonstrate, is ideally suited to the typical, finite-temperature 'impurity' problems that describe open quantum system connected to multiple environments.
Out-of-equilibrium dynamics of non-integrable Hamiltonian many-body quantum systems are characterized by highly entangled wave functions. Near-maximal entanglement arises in systems exhibiting thermalization or pre-thermalization, where the system converges to a steady state with a fixed energy density. Classical simulation of the time dependence of such wave functions requires exponential resources. However, typical computations aim to estimate expectation values of local operators and correlation functions to some expected precision. For thermalizing systems at sufficiently high energy densities, such computations can be done without storing the full wave function by instead simulating the evolution of the local operator, which requires significantly fewer resources. Nonetheless, constructing such resource-efficient classical algorithms remains a challenge for intermediate energy densities, where simulating both the wave function and operator evolution is costly. In this paper, we propose a heuristic approach to accelerate the convergence of Matrix Product State (MPS) simulations of expectation values, applicable across a broad range of energy densities. We estimate the desired observables by rescaling the MPS results at low bond dimensions with a factor that depends on the fidelity of the MPS wave function. Using this technique, we simulated the dynamics of the two-dimensional Transverse-Field Ising Model (TFIM) on a 7×8 grid with periodic boundary conditions, using a maximum bond dimension of χ=4096 on a single A100 GPU, as well as the dynamics of the two-dimensional XY model on grids of size up to 9×9. We compare our TFIM results to similar simulations on a digital quantum processor [R. Haghshenas et al., Nature 653, 56 (2026)], demonstrating excellent agreement and confirming the predictive power of our method.
We present a unified and highly numerically efficient formalism for the simulation of quantum dynamics of complex molecular systems, which takes into account both temperature effects and static disorder. The methodology is based on the thermo-field dynamics formalism, and Gaussian static disorder is included into simulations via auxiliary bosonic operators. This approach, combined with the tensor-train/matrix-product state representation of the thermalized stochastic wave function, is applied to study the effect of dynamic and static disorders in charge-transfer processes in model organic semiconductor chains employing the Su-Schrieffer-Heeger (Holstein-Peierls) model Hamiltonian.
Modeling open quantum systems---quantum systems coupled to a bath---is of value in condensed-matter theory, cavity quantum electrodynamics, nanosciences, and biophysics. The real-time simulation of open quantum systems was advanced significantly by the recent development of chain mapping techniques and the use of matrix product states that exploit the intrinsic entanglement structure in open quantum systems. The computational cost of simulating open quantum systems, however, remains high when the bath is excited to high-lying quantum states. We develop an approach to reduce the computational costs in such cases. The interaction representation for the open quantum system is used to distribute excitations among the bath degrees of freedom so that the occupation of each bath oscillator is ensured to be low. The interaction picture also causes the matrix dimensions to be much smaller in a matrix product state of a chain-mapped open quantum system than in the Schr\"odinger picture. Using the interaction representation accelerates the calculations by one to two orders of magnitude over the existing matrix-product-state method. In the regime of strong system-bath coupling and high temperatures, the speedup can be as large as three orders of magnitude. The approach developed here is especially promising to simulate the dynamics of open quantum systems in high-temperature and strong-coupling regimes.
Matrix-product states (MPS) have proven to be a versatile ansatz for modeling quantum many-body physics. For many applications, and particularly in one-dimension, they capture relevant quantum correlations in many-body wave functions while remaining tractable (polynomial scaling) to store and manipulate on a classical computer. This has motivated researchers to also apply the MPS ansatz to machine learning (ML) problems where capturing complex correlations in datasets is also a key requirement. Here, for the first time, we develop an MPS-based algorithm, MPSTime, for learning a joint probability distribution underlying a time-series dataset, and show how it can be used to tackle important time-series ML problems, including classification and imputation. MPSTime can efficiently learn complicated time-series probability distributions directly from data, requires only moderate maximum MPS bond dimension χ max (with values for our applications ranging between χ max = 20 − 160 ), and can be trained for both classification and imputation tasks under a single logarithmic loss function. Using synthetic and publicly available real-world datasets—spanning applications in medicine, energy, and astronomy—we demonstrate performance competitive with state-of-the-art ML approaches, but with the key advantage of encoding the full joint probability distribution learned from the data, which is useful for analyzing and interpreting its underlying structure. This manuscript is supplemented with the release of a publicly available code package that implements our approach, and can be used to reproduce all presented results. The effectiveness of the MPS-based ansatz for capturing complex correlation structures in time-series data makes it a powerful foundation for tackling challenging time-series analysis problems across science, industry, and medicine.
Efficient simulation of open quantum systems is central to modeling noisy quantum hardware and many-body dynamics. In trajectory-based tensor network methods, cost is often associated with trajectory-level quantities such as entanglement growth or bond dimension. However, the total cost of a fixed-accuracy simulation also depends on statistical sampling, and the interplay between per-trajectory complexity and sampling effort remains poorly understood. Here we introduce a cost-resolved framework for matrix product state (MPS)-based quantum trajectory simulations that decomposes total cost into memory per trajectory, runtime per trajectory, and sampling effort. We show that physically equivalent stochastic unravelings of the same Lindblad dynamics do not necessarily reduce total cost, but instead redistribute cost between trajectory complexity and statistical convergence. This trade-off is quantified by two dimensionless inflation factors: a bond dimension inflation α and a sampling inflation κ, which together determine the preferred unraveling under hardware-dependent memory and parallelism constraints. We provide a practical protocol for extracting (α,κ) from modest pilot simulations and demonstrate it using benchmarks across multiple noise channels. The resulting decision maps show that the computationally favorable unraveling can change with noise strength, time-step resolution, system size, and available parallelism. These results establish unraveling choice as a hardware-aware simulation design problem rather than an intrinsic optimization of trajectory entanglement alone.