Determining the steady state of an open quantum system is crucial for characterizing quantum devices and studying various physical phenomena. Often, computing a single steady state is insufficient, and it is necessary to explore its dependence on multiple external parameters. In such cases, calculating the steady state independently for each combination of parameters quickly becomes intractable. Perturbation theory (PT) can mitigate this challenge by expanding steady states around reference parameters, minimizing redundant computations across neighboring parameter values. However, PT has two significant limitations: it relies on the pseudo-inverse – a numerically costly operation – and has a limited radius of convergence. In this work, we remove both of these roadblocks. First, we introduce a variational perturbation theory (VPT) and its multipoint generalization that significantly extends the radius of convergence even in the presence of non-analytic effects such as dissipative phase transitions. Then, we develop two numerical strategies that eliminate the need to compute pseudo-inverses. The first relies on a single LU decomposition to efficiently construct the steady state within the convergence region, while the second reformulates VPT as a Krylov space recycling problem and uses preconditioned iterative methods. We benchmark these approaches across various models, demonstrating their broad applicability and significant improvements over standard PT.
Determining the steady state of an open quantum system is crucial for characterizing quantum devices and studying various physical phenomena. Often, computing a single steady state is insufficient, and it is necessary to explore its dependence on multiple external parameters. In such cases, calculating the steady state independently for each combination of parameters quickly becomes intractable. Perturbation theory (PT) can mitigate this challenge by expanding steady states around reference parameters, minimizing redundant computations across neighboring parameter values. However, PT has two significant limitations: it relies on the pseudo-inverse – a numerically costly operation – and has a limited radius of convergence. In this work, we remove both of these roadblocks. First, we introduce a variational perturbation theory (VPT) and its multipoint generalization that significantly extends the radius of convergence even in the presence of non-analytic effects such as dissipative phase transitions. Then, we develop two numerical strategies that eliminate the need to compute pseudo-inverses. The first relies on a single LU decomposition to efficiently construct the steady state within the convergence region, while the second reformulates VPT as a Krylov space recycling problem and uses preconditioned iterative methods. We benchmark these approaches across various models, demonstrating their broad applicability and significant improvements over standard PT.
We present a general nonperturbative method to determine the exact steady state of open quantum systems under perturbation. The method works for systems with a unique steady state and the perturbation may be time-independent or periodic, and of arbitrarily large amplitude. Using the Drazin inverse and a single diagonalization, we construct an operator that generates the entire dependence of the steady state on the perturbation parameter. The approach also enables exact analytic operations—such as differentiation, integration, and ensemble averaging—with respect to the parameter, even when the steady state is computed numerically. We apply the method to three nontrivial open quantum systems, showing that it achieves exact results, with a computational speedup of one to several orders of magnitude for calculations requiring large sampling, compared to previous approaches.
Abstract In open quantum systems, the Liouvillian gap characterizes the relaxation time toward the steady state. However, accurately computing this quantity is notoriously difficult due to the exponential growth of the Hilbert space and the non-Hermitian nature of the Liouvillian superoperator. In this work, we propose a variational quantum algorithm for efficiently estimating the Liouvillian gap. By utilizing the Choi–Jamiołkowski isomorphism, we reformulate the problem as finding the first excitation energy of an effective non-Hermitian Hamiltonian. Our method employs variance minimization with an orthogonality constraint to locate the first excited state and adopts a two-stage optimization scheme to enhance convergence. Moreover, to address scenarios with degenerate steady states, we introduce an iterative energy-offset scanning technique. Numerical simulations on the dissipative XXZ model confirm the accuracy and robustness of our algorithm across a range of system sizes and dissipation strengths. These results demonstrate the promise of variational quantum algorithms for simulating open quantum many-body systems on near-term quantum hardware.
We discuss the numerical solution methods available when solving for the steady-state density matrix of a time-independent open quantum optical system, where the system operators are expressed in a suitable basis representation as sparse matrices. In particular, we focus on the difficulties posed by the non-Hermitian structure of the Lindblad super operator, and the numerical techniques designed to mitigate these pitfalls. In addition, we introduce a doubly iterative inverse-power method that can give reduced memory and runtime requirements in situations where other iterative methods are limited due to poor bandwidth and profile reduction. The relevant methods are demonstrated on several prototypical quantum optical systems where it is found that iterative methods based on iLU factorization using reverse Cuthill-Mckee ordering tend to outperform other solution techniques in terms of both memory consumption and runtime as the size of the underlying Hilbert space increases. For eigenvalue solving, Krylov iterations using the stabilized bi-conjugate gradient method outperform generalized minimal residual methods. In contrast, minimal residual methods work best for solvers based on direct LU decomposition. This work serves as a guide for solving the steady-state density matrix of an arbitrary quantum optical system, and points to several avenues of future research that will extend the applicability of these classical algorithms in absence of a quantum computer.
The theory of open quantum systems lays the foundations of a substantial part of modern research in quantum science and engineering. Rooted in the dimensionality of their extended Hilbert spaces, the high computational complexity of simulating open quantum systems calls for the development of strategies to approximate their dynamics. In this paper, we present an approach for tackling open quantum system dynamics. We simulate the dynamics of the Liouvillian superoperator using a forward-backward trapezoid method and find the steady-state via a variational formulation. We make use of a probabilistic formulation of quantum physics based on a positive operator-valued measure (POVM) in combination with neural autoregressive models, which bring significant algorithmic flexibility due to their efficient sampling and tractable density. We introduce improved ansatzs, String States, which partially restore the symmetry of the neural auto-regressive model and improve the description of local correlations. We benchmark our approaches on prototypical one and two-dimensional systems, finding results which closely track the exact solution and achieve higher accuracy in comparison to a recently proposed method based on restricted Boltzmann machines. We anticipate this approach will be widely applicable to evolving density matrices in various contexts.
Open quantum systems have drawn attention over decades due to its applicability in the foundation of theoretical physics, e.g. statistical mechanics, quantum mechanics and condensed matter physics. The dynamics of open quantum systems has been described as separate entities from their surrounding environment that consists of a very large number of modes, somehow coupled to the mode of the system. Even though the exact solution of the dynamical behavior of the system is impossible to calculate, we obtain a tentative solution using the crucial Markov approximation. The input-output formalism of the quantum Langevin equation (QLE) has been considered as a useful tool which provides a semi-classical description of the dynamics of the system, whereas the master equation provides a complete picture of the dynamics of the system expressed in terms of density matrix. While studying the dynamics of nonlinear system/environment coupling using QLE, we see, for a small value of external field, that the steady state system field does not change much from the steady state field obtained in the absence of nonlinear dissipation. However, in a case where a stronger external field is applied, we see that the deviation becomes substantial from the solution of linear system. We also see that the nonlinear coupling introduces significant difference in the cavity fluctuation spectrum. The description, therefore, provides a potential explanation of parametric effects in terms of nonlinear dissipation phenomena associated with the nonlinear coupling. Even though the theories developed in the context of open quantum systems have proven to be powerful tools, they do not provide a satisfactory platform to be implemented on non-linear Hamiltonians. We often approximate it by linearizing over nonlinear steady state field amplitude, and therefore, the interesting effects are often overlooked. The limitation of the analytics provokes us to simulate open quantum dynamics numerically. The numerical method consists of transformation of the environmental degrees of freedom to a one-dimensional many-body chain, and the computational technique includes numerical diagonalization and renormalization process. The time-adaptive density matrix renormalisation group (t-DMRG) is known as one of the most powerful techniques for the simulation of strongly-correlated many-body quantum systems. In this thesis, along with the theoretical modeling, we implement DMRG numerical scheme for the simulation of canonical S/B model by mapping it to one-dimensional harmonic chain with nearest neighbor interactions, and use the method to investigate the dynamics of the free dissipative system. The thermalization of open quantum systems is also studied by generating minimally entangled typical thermal states (METTS) through imaginary time evolution, and real-time evolving an empty system in the presence of the thermal bath. Further, we simulate coherently driven free dissipative Kerr nonlinear system numerically using Euler's method by solving Heisenberg equation of motion and t-DMRG algorithm, and demonstrate how the numerical results are analogous to classical bistability. By comparing with analytics, we see that the DMRG numerics is analogous to the quantum-mechanical exact solution obtained by mapping the equation of motion of the density matrix of the system to a Fokker-Plank equation. The comparison between two different numerical techniques shows that the semi-classical Euler's method determines the dynamics of the system field of one among two coherent branches, whereas DMRG numerics gives the superposition of both of them. Hence, DMRG-determined time dynamics undergoes generating non-classical states. Our approach of dealing with nonlinearity represents an important contribution in the developments of technique to study the dynamical and steady-state behavior of open quantum systems, which is a fundamental aspect of quantum physics.
Non-equilibrium quantum phase transitions (NEQPTs) represent a frontier in condensed matter physics and quantum statistical mechanics, characterized by abrupt changes in the steady-state properties of open quantum systems or by critical dynamics in closed systems far from thermal equilibrium. Unlike their equilibrium counterparts, NEQPTs are profoundly influenced by driving forces, dissipation, and the specific dynamics of quantum quenches. Understanding these complex phenomena requires theoretical frameworks that can accurately capture both strong correlations and real-time evolution in the presence of environmental interactions. Traditional theoretical tools, such as purely variational methods or purely perturbative expansions, often fall short when dealing with the intricate interplay of coherent and incoherent processes characteristic of NEQPTs. Variational approaches, while adept at capturing ground state properties and short-time dynamics within a chosen ansatz, struggle with long-time thermalization and the inclusion of subtle quantum and thermal fluctuations beyond the variational manifold. Conversely, perturbative methods, particularly those based on the Keldysh formalism, can systematically include interactions and environmental effects but typically rely on a small expansion parameter, limiting their applicability to weakly interacting or weakly driven regimes. This paper proposes a novel "Variational-Perturbative Unification" framework designed to overcome these limitations. The core idea is to combine the strengths of both approaches: using a variational ansatz to capture the dominant non-equilibrium many-body state or dynamics, and then systematically incorporating quantum and thermal fluctuations, as well as dissipation, through a perturbative expansion around this optimized variational state. This hybrid methodology aims to provide a robust and accurate description of NEQPTs, offering insights into critical exponents, universality classes, and the dynamic pathways to steady states. We explore the theoretical underpinnings of this unified approach, discuss its potential applications to driven-dissipative systems and quantum quenches, and highlight how it can reveal novel non-equilibrium phenomena.