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

TL;DR

The paper addresses the costly problem of computing steady states for open quantum systems under parameter variation by introducing variational perturbation theory (VPT), which recasts perturbative corrections as a low-dimensional variational problem to extend the convergence radius and avoid the Moore–Penrose inverse. It proposes two model-independent strategies—LU-based reuse and preconditioned Krylov methods—to compute steady states and their gradients efficiently, and extends PT to multipoint VPT (m-VPT) to handle non-analytic behavior near dissipative phase transitions. Through benchmarks on driven Kerr resonators, dissipative cat, and dissipative XYZ models, the authors demonstrate substantial speedups and robust phase-diagram and parameter-estimation capabilities. The framework is model-agnostic, compatible with other reduction techniques, and poised for extensions to time-domain simulations and more complex systems. Overall, VPT provides a practical, scalable approach to exploring open quantum systems across parameter spaces with improved convergence and computational efficiency.

Abstract

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.

Figures (6)

• Figure 1: Depiction of the working principle of standard perturbation theory (PT), variational PT (VPT), and of the perturbative recursion relation. We are interested in finding steady states of a parameterized Liouvillian $\mathcal{L}(\varepsilon) = \mathcal{L}_0 + \varepsilon \mathcal{L}_1$ for many values of $\varepsilon$ [or, more generally, $\mathcal{L}(\varepsilon_1, \varepsilon_2, \dots \varepsilon_n) = \mathcal{L}_0 + \sum \varepsilon_j \mathcal{L}_{j}$]. (a) Sketch of the parameter space. At the point $\varepsilon =0$, the steady state $\hat{\rho}_{\rm ss}(\varepsilon=0)$ of $\mathcal{L}(\varepsilon=0) = \mathcal{L}_0$ is computed. To compute the steady state in neighbouring points in the parameter space, we resort to PT and find the set of matrices $\hat{\rho}_{\rm ss}^{(n)}$ defined in Eq. (\ref{['Eq:recursion_relation']}). We develop methods that do that efficiently for (b) 1D parameter spaces, (c) 2D parameter spaces and, more generally, for arbitrarily high dimensional parameter space. These methods are based on either re-applying the $LU$ decomposition used to compute $\hat{\rho}_{\rm ss}(\varepsilon=0)$, or using approximate methods to build the basis of matrices defining the steady state. (d) Having obtained the series of perturbation matrices $\hat{\rho}_{\rm ss}^{(n)}$, one can construct the state $\hat{\rho}_{\rm ss}(\varepsilon)$ by weighting each matrix according to the power series in Eq. (\ref{['Eq:PT_coeff']}) within the radius of convergence. (e) Allowing for a more expressive ansatz in the form of Eq. (\ref{['Eq:PT_coeff_general']}), one can use $\hat{\rho}_{\rm ss}^{(n)}$ to describe the solution on a much wider range of parameters. This makes VPT more efficient than PT, as the "heavy" operation of computing the steady state at one point via exact numerical methods becomes less frequent.
• Figure 2: Comparison of standard PT and VPT in the study of the driven-dissipative Kerr resonator described by Eq. (\ref{['Eq:Kerr_resonator']}) and sketched in (a). (b) Average photon number $\expval{\hat{a}^\dagger \hat{a}}$ as a function of the detuning $\Delta$ determined using an exact solution (black line), standard PT (yellow line), and VPT (blue line). Both PTs have been performed starting from $\Delta =0$, as indicated by the dashed vertical line. The solid background denotes the region where the perturbative expansion has converged up to $\| \mathcal{L}(\varepsilon) \hat{\rho}_{\rm ss}^{\rm PT}(\varepsilon) \|<10^{-2}$. In the yellow region, both standard PT and VPT are convergent. In the blue region, only VPT reached convergence. (c) Error $\| \mathcal{L}(\varepsilon) \hat{\rho}_{\rm ss}^{\rm PT}(\varepsilon) \|$ as a function of the maximal order of perturbation $M$ for standard PT. (d) Same as (b), but for VPT. (e-h) Coefficients obtained by standard [c.f. Eq. (\ref{['Eq:PT_coeff']})] and variational [c.f. Eq. (\ref{['Eq:PT_coeff_general']})] PTs. Parameters: $K/\kappa = 10$ and $\varepsilon/\kappa = 10$. Fock space truncation: 30 photons.
• Figure 3: Comparison of the performance of perturbation theory methods to compute the steady of the driven-dissipative Kerr resonator defined in Eq. (\ref{['Eq:Kerr_resonator']}) in a two-dimensional parameter space. In (a) we plot the average photon number $\langle a^\dagger a\rangle$ in the steady state as a function of detuning $\Delta$ and drive strength $F$ along with the boundaries of the convergence regions of multipoint VPT. In (b) we plot the points where we exactly computed the steady state through LU decomposition. At each point we computed a grid of perturbative corrections up to $M = (10, 10)$, totaling 121 perturbation vectors per point, and set the error tolerance to $\| \mathcal{L}(\varepsilon) \hat{\rho}_{\rm ss}^{\rm PT}(\varepsilon) \|<10^{-7}$. Parameters: $K/\kappa = \frac{1}{2}$. Fock space truncation: 70 photons.
• Figure 4: Performance of perturbation theory methods with increasing perturbation order $M$. We study the same system as in Fig. 3 and apply (V)PT at 40 randomly selected points in parameter space. For each resulting convergence region, we compute the optimal low-rank basis (see Appendix \ref{['Appendix:SVD']}. In both panels, blue violins correspond to VPT and yellow violins to PT, with scatter points indicating individual samples. (a) Violin plot of the convergence region area as a function of the VPT basis size $(M+1)^2$. (b) Violin plot of $(M+1)^2 / M_\mathrm{SVD}$ as a function of $(M+1)^2$, where $(M+1)^2$ is the number of perturbation vectors used in (V)PT, and $M_\mathrm{SVD}$ is the size of the optimal low-rank basis. A higher ratio indicates a more efficient representation relative to the optimal basis. Violin plots illustrate data distributions, with the width representing the density of data points at each value, estimated using a Gaussian kernel.
• Figure 5: Parameter estimation using VPT and gradient-based optimization for a memory-buffer system sketched in (a) and described by the Liouvillian in Eq. (\ref{['Eq:Liouvillian_cat']}). (b) Synthetic data of the steady state reflection coefficient $S_{12}$ as a function of $\Delta_a$ and $\Delta_b$. To mimic experimental imperfections, we add Gaussian noise with average $\mu =0$ and standard deviation $\sigma = 0.02$. (c-e) Reflection coefficient at various stages of the optimization. The algorithm starts from a random guess of $g_2$, $F$, and $K_a$ and uses L-BFGS algorithm to optimize them, assuming the remaining parameters are known. In (f) we plot the loss as a function of the iteration number, while (g-i) show how the parameters evolve along the optimization process. The dashed lines indicate the true parameters used to generate the data. Parameters (in MHz): $g_2 = 2$, $F = 2$, $K_a = 0.1$, $K_b = 0.3$, $\kappa_a = 0.1$, $\kappa_b = 10$. We fix the Fock space truncation at $10$ photons for mode $\hat{a}$ and $6$ for $\hat{b}$.