Variational Perturbation Theory in Open Quantum Systems for Efficient Steady State Computation

Determining the steady state of an open quantum system is essential for characterizing quantum devices and understanding diverse physical phenomena. In many cases, however, computing a single steady state is insufficient—one must explore its dependence on several external parameters. Computing these steady states independently for every parameter combination quickly becomes intractable. Perturbation theory alleviates this challenge by expanding steady states around reference parameters, thereby avoiding redundant computations across neighboring points. Yet, it faces two major obstacles: it depends on the pseudo-inverse, a numerically expensive operation, and it has a limited radius of convergence. In this work, we overcome both limitations. We first introduce variational perturbation theory and its multipoint generalization, which greatly extend the convergence radius even in the presence of non-analytic behavior such as dissipative phase transitions. Next, we develop two complementary numerical strategies that eliminate the need to compute pseudo-inverses. The first constructs steady states efficiently within the convergence region using a single LU decomposition, while the second reformulates VPT as a Krylov subspace recycling problem solved via preconditioned iterative methods. We benchmark these approaches across several models, demonstrating their wide applicability and substantial performance gains over standard PT.