A Variance Reduction Technique for the Stochastic Liouville-von Neuman Equation

TL;DR

This paper introduces a variance reduction sampling strategy for the stochastic Liouville-von Neumann equation, significantly improving simulation efficiency for quantum systems with dissipative interactions.

Contribution

It proposes a novel sampling method that reduces non-unitary terms via convex optimization, enhancing scalability with dissipative interaction strength.

Findings

Improved scaling with dissipation strength

Reduced variance in stochastic simulations

Enhanced efficiency over previous methods

Abstract

The Stochastic Liouville-von Neumann equation provides an exact numerical simulation strategy for quantum systems interacting with Gaussian reservoirs [J.T. Stockburger & H. Grabert, PRL 88, 170407 (2002)]. Its scaling with the extension of the time interval covered has recently improved dramatically through time-domain projection techniques [J.T. Stockburger, EPL 115, 40010 (2016)]. Here we present a sampling strategy which results in a significantly improved scaling with the strength of the dissipative interaction, based on reducing the non-unitary terms in sample propagation through convex optimization techniques.

Figures

Tables

Table 1. Table 1: Overview of the correspondence between Gaussian quantum fluctuations of an environment and stochastic processes with matched statistics.

	trace reduction	stochastic reduction
state before reduction:	$\rho$ in global Liouville space	random $\rho$ in system space
reservoir “forces”:	$B_{+}$, $B_{-}$ superoperators	$\xi \left(t\right)$, $\nu \left(t\right)$ random $\in ℂ$
reduction operation:	partial trace ${\mathrm{tr}}_R$	expectation value
	with weight ${\rho }_R$	w.r.t. noise probability density
	$⟨B_{+}\left(t\right)B_{+}\left(t^{\prime }\right)⟩$	$⟨\xi \left(t\right)\xi \left(t^{\prime }\right)⟩$
	$⟨B_{+}\left(t\right)B_{-}\left(t^{\prime }\right)⟩$	$⟨\xi \left(t\right)\nu \left(t^{\prime }\right)⟩$
	$⟨B_{-}\left(t\right)B_{+}\left(t^{\prime }\right)⟩=0$	$⟨\nu \left(t\right)\xi \left(t^{\prime }\right)⟩=0$
	$⟨B_{-}\left(t\right)B_{-}\left(t^{\prime }\right)⟩=0$	$⟨\nu \left(t\right)\nu \left(t^{\prime }\right)⟩=0$
	time ordering $t>t^{\prime }$	time ordering “by hand”

Table 2. Table 2: Performance characteristics without (left of double line) and with noise optimization. The number of samples  nsampnote results from demanding an absolute error of 0.01 0.01 0.01 for the simulation result ⟨ σ x ⟩ delimited-⟨⟩ subscript 𝜎 𝑥 \langle\sigma_{x}\rangle . All data are for a symmetric spin-boson system with K = 0.05 𝐾 0.05 K=0.05 .

	$\beta \Delta =5$	$\beta \Delta =5$	$\beta \Delta =2$	$\beta \Delta =1$	$\beta \Delta =0.5$
variance growth rate ($\Delta$)	$0.20$	$0.10$	$8.6\cdot {10}^{-2}$	$5.5\cdot {10}^{-2}$	$2.5\cdot {10}^{-2}$
variance at $\Delta t=40$	$3\cdot {10}^3$	$5.2\cdot {10}^2$	$1.8\cdot {10}^1$	$5.5\cdot {10}^0$	$1.7\cdot {10}^0$
extrapolated number
of samples at $\Delta t=40$	$3\cdot {10}^7$	$5.2\cdot {10}^6$	$1.8\cdot {10}^5$	$5.5\cdot {10}^4$	$1.7\cdot {10}^4$

Equations

$$\dot{\rho }(t)=i{A}_{+}(t)B_{-}(t)\rho (t)+i{A}_{-}(t)B_{+}(t)\rho (t).$$

$$\rho (t)={\exp }_{>}\left(i{\int }_0^{t}ds(A_{+}(s)\otimes B_{-}(s)+A_{-}(s)\otimes B_{+}(s))\right)\rho (0).$$

$${\rho }_S(t)$$

$$\displaystyle=\exp_{>}\bigg{(}i\int_{0}^{t}ds\int_{0}^{s}ds^{\prime}(A_{-}(s)A_{-}(s^{\prime})\big{\langle}B_{+}(s)B_{+}(s^{\prime})\big{\rangle}$$

$$\displaystyle\qquad\qquad\qquad\qquad\qquad+A_{-}(s)A_{+}(s^{\prime})\big{\langle}B_{+}(s)B_{-}(s^{\prime})\big{\rangle})\bigg{)}\rho_{S}(0).$$

$$⟨\xi (t)\xi (t^{\prime })⟩$$

$$⟨\xi (t)\nu (t^{\prime })⟩$$

$$⟨\nu (t)\nu (t^{\prime })⟩$$

$$\dot{\tilde{\rho }}={\mathcal{L}}_S\tilde{\rho }+i\xi (t)A_{-}\tilde{\rho }+i\nu (t)A_{+}\tilde{\rho },$$

$$⟨{\xi }^{(l)}(t){\xi }^{(l)}(t^{\prime })⟩$$

$$⟨{\xi }^{(s)}(t)\nu (t^{\prime })⟩$$

$$\mathsf{\Sigma}(\omega)=F\left[\left(\begin{array}[]{cccc}\langle{\xi^{\prime}}(t){\xi^{\prime}}(t^{\prime})\rangle&\langle{\xi^{\prime}}(t){\xi^{\prime\prime}}(t^{\prime})\rangle&\langle{\xi^{\prime}}(t){\nu^{\prime}}(t^{\prime})\rangle&\langle{\xi^{\prime}}(t){\nu^{\prime\prime}}(t^{\prime})\rangle\\ \langle{\xi^{\prime\prime}}(t){\xi^{\prime}}(t^{\prime})\rangle&\langle{\xi^{\prime\prime}}(t){\xi^{\prime\prime}}(t^{\prime})\rangle&\langle{\xi^{\prime\prime}}(t){\nu^{\prime}}(t^{\prime})\rangle&\langle{\xi^{\prime\prime}}(t){\nu^{\prime\prime}}(t^{\prime})\rangle\\ \langle{\nu^{\prime}}(t){\xi^{\prime}}(t^{\prime})\rangle&\langle{\nu^{\prime}}(t){\xi^{\prime\prime}}(t^{\prime})\rangle&\langle{\nu^{\prime}}(t){\nu^{\prime}}(t^{\prime})\rangle&\langle{\nu^{\prime}}(t){\nu^{\prime\prime}}(t^{\prime})\rangle\\ \langle{\nu^{\prime\prime}}(t){\xi^{\prime}}(t^{\prime})\rangle&\langle{\nu^{\prime\prime}}(t){\xi^{\prime\prime}}(t^{\prime})\rangle&\langle{\nu^{\prime\prime}}(t){\nu^{\prime}}(t^{\prime})\rangle&\langle{\nu^{\prime\prime}}(t){\nu^{\prime\prime}}(t^{\prime})\rangle\end{array}\right)\right]$$

$${\Sigma }_{{\xi }^{\prime }{\xi }^{\prime }}+i{\Sigma }_{{\xi }^{\prime }{\xi }^{\prime \prime }}+i{\Sigma }_{{\xi }^{\prime \prime }{\xi }^{\prime }}-{\Sigma }_{{\xi }^{\prime \prime }{\xi }^{\prime \prime }}$$

$${\Sigma }_{{\xi }^{\prime }{\nu }^{\prime \prime }}+{\Sigma }_{{\xi }^{\prime \prime }{\nu }^{\prime }}$$

$${\Sigma }_{{\xi }^{\prime }{\nu }^{\prime }}+{\Sigma }_{{\xi }^{\prime \prime }{\nu }^{\prime \prime }}$$

$${\Sigma }_{{\nu }^{\prime }{\nu }^{\prime }}+i{\Sigma }_{{\nu }^{\prime }{\nu }^{\prime \prime }}+i{\Sigma }_{{\nu }^{\prime \prime }{\nu }^{\prime }}-{\Sigma }_{{\nu }^{\prime \prime }{\nu }^{\prime \prime }}$$

$$\displaystyle\left(\begin{array}[]{ll}\Sigma_{\xi^{\prime}\xi^{\prime}}&\Sigma_{\xi^{\prime}\nu^{\prime\prime}}\\ \Sigma_{\nu^{\prime\prime}\xi^{\prime}}&\Sigma_{\nu^{\prime\prime}\nu^{\prime\prime}}\end{array}\right)$$

$$A$$

$$C$$

Full text

1 Introduction

The notion of an open quantum system originates from the embedding of a system of interest into a larger environment, where system and environment together are governed by ordinary unitary quantum evolution. Different abstractions and techniques, both formal and computational have been put forward to formulate a reduced dynamics of the system of interest without explicit reference to the environment beyond its initial preparation.

Completely positive trace-preserving channels provide a very broadly applicable formal description of the time evolution of open quantum systems alick87 . The proof of their existence is constructive; it does not, however, in itself provide a dynamical law governing the reduced description. For a quantum system with discrete energy spectrum, the respective quantum master equation can be derived within the constraints of a combined Born-Markov-secular approximation davie74 ; breue02 .

The dissipator resulting from this depends not only on reservoir and coupling properties, but also on the energy spectrum and the eigenstates of the system. For complicated spectra, its determination can be cumbersome. The simpler alternative of adapting separate dissipators to individual parts of a system (“local” Lindblad operators) may result in dynamics which contradicts elementary properties of thermal environments levy14 ; stock17 , e.g., unphysical heat flow against a thermal gradient. Similar problems arise in the context of driven quantum systems alick06 ; schmi11 .

The formal description of environmental effect through influence functionals feynm63 ; weiss08b has the attractive feature of describing a dissipation mechanism without direct reference to any system properties beyond specifying the coupling Hamiltonian. It is fully applicable in the context of any driving. In fact, it is completely “agnostic” to the nature of the free system dynamics, since it keeps the given separable structure of the interaction Lagrangian. Direct Monte Carlo evaluation of real-time path integrals is feasible egger00 ; muhlb04 but numerically expensive in regimes without rapid dephasing by the environment (dynamical sign problem). Equations of motion equivalent to the influence functional approach can be found only at the price of introducing a large hierarchy of auxiliary states tanim91 ; tanim14 , since the weights assigned to paths by influence functionals are non-local in time.

On the other hand, the stochastic unraveling of the influence functional can restore this locality, at least on the level of individual realizations of auxiliary trajectories (which can be interpreted as Gaussian colored noise): The Feynman-Vernon influence functional feynm63 has essentially the same formal structure as the generating functional of Gaussian noise, hence it can be viewed as the noise average over exponentials of a time-local action functional stock02 . This idea has led to a variety of stochastic propagation schemes cao96 ; strun96 ; stock02 ; shao04 ; tanim06 . The corresponding numerical methods use independently drawn samples of the auxiliary Gaussian noise trajectories. For each noise sample, the time evolution of the reduced density matrix is obtained by integration of a Liouville equation with random terms, in a similar manner as a direct simulation of a generalized Langevin equation. The physical density matrix of the system is recovered by taking the average of these sample density matrices in the limit of a large number of samples. It is to be noted that the sample density matrices are not physical quantum states; only their expectation value with respect to the probability measure of the Gaussian noise is to be identified with the reduced density matrix of the system.

The complete unraveling of Gaussian influence functionals introduced by Stockburger and Grabert stock02 is valid for arbitrary spectral properties and temperatures of the reservoir. Its direct practical application is sometimes hampered by rapid growth of the sample variance when long propagation intervals are covered. This problem has recently been solved for the near-universal case of a finite (not necessarily small) correlation time of the free reservoir fluctuations stock16a . The required number of samples has thus been reduced – by orders of magnitude in some parameter regimes. For combinations of moderate to strong coupling and long reservoir memory times, the number of samples may be sufficiently high to require the use of parallel computational resources.

The present paper addresses this remaining problem by optimization of the noise correlation functions, which are not entirely determined by the unraveling procedure. Section 2 gives a brief overview of the unraveling of propagating functions based on influence functionals and their transformation into a stochastic Liouville-von Neumann (SLN) equation, including a discussion of issues related to sampling statistics. In section 3 a strategy is developed which improves sampling statistics through minimizing the power of problematic noise components, and numerical examples are provided. Section 4 comprises conclusions and an outlook.

2 Stochastic unraveling and SLN equation

The interacting dynamics of a system S coupled to an environment (reservoir) R is governed by a Hamiltonian which can be partitioned as $H=H_{S}+H_{R}+H_{I}$, where the indices S and R stand for terms which act exclusively on system and reservoir degrees of freedom. For the interaction we assume separability, $H_{I}=-A\otimes B$, where $A$ and $B$ act on system and reservoir.

We choose an initially thermal state of the environment (uncorrelated with the system) with the intent of tracing out the reservoir degrees of freedom from the correlated state arising from the dynamics. The unitary propagation of the global system-plus-reservoir density matrix is the appropriate conceptual starting point for this approach. Instead of using path integrals for this purpose stock02 , equivalent time-ordered exponentials will be considered here. These can be considered generating functionals of non-commuting random variables, and they simplify to Gaussian expressions in exactly the same cases for which one obtains Gaussian influence functionals kubo62 .

A most compact form of the dynamics hbarone is obtained in the interaction representation of the Liouville dynamics,

[TABLE]

For reasons which will become clear later, a choice was made to decompose the interaction Liouvillian into superoperators for reservoir $B_{\pm}$, and system $A_{\pm}$ defined through (anti-)commutators, $B_{-}=[B,\cdot]$ and $B_{+}=\frac{1}{2}\{B,\cdot\}$ etc.

The formal solution of this dynamics is the time-ordered exponential

[TABLE]

We now consider the case of Gaussian statistics of $B$ and an initially factorizing state, $\rho_{S}(0)\otimes\rho_{R}(0)$. Using the notation $\langle\cdot\rangle=\mathop{\mathrm{tr}_{\mathrm{R}}}\left\{\cdot\rho_{R}(0)\right\}$ for the partial-trace averaging procedure, the reduced density matrix is now formally given by

[TABLE]

Of the four terms formally arising in a bivariate Gaussian characteristic function, two are identically zero here because they are traces over a commutator. A key observation at this point is the following: Had we started with a different averaging procedure – replacing $B_{+}$ and $B_{-}$ by $c$-number Gaussian noise – we would have arrived at a similar expression.

The conditions

[TABLE]

are necessary and sufficient for the exact identification of the open-system dynamics (3) and stochastic propagation with the substitutions $B_{+}\to\xi$ and $B_{-}\to\nu$. On the left hand side, the angle brackets stand for expectation values with respect to the noise statistics. In practice, the stochastic equivalent of eq. (3) is evaluated by returning to the Schrödinger picture and averaging solutions of the SLN equation

[TABLE]

where $\tilde{\rho}$ is the sample density matrix associated with any particular realization of $\xi(t)$ and $\nu(t)$.

Averaging $\tilde{\rho}$ over samples provides an estimate of the physical density matrix $\rho_{S}$. Since noise samples are generated independently, an accurate estimate of the statistical error of this approach can be determined. In this context, we also note that parallelization of the simulation method is trivial. An overview of the quantum-stochastic correspondence is given in Table 1.

It is now crucial to observe that eqs. (4–6) cannot be obeyed by real-valued noise. Only when extending at least $\nu$ to values in the complex plane can these equations hold: $\xi$ and $\nu$ are random variables in the sense of ordinary Kolmogorov probability. When extending both $\xi=\xi^{\prime}+i\xi^{\prime\prime}$ and $\nu=\nu^{\prime}+i\nu^{\prime\prime}$ to complex values, the correlations of their four real components are not completely determined by eqs. (4–6). Additional non-physical correlation functions, in complex notation, $\langle\xi^{*}(t)\xi(t^{\prime})\rangle$, $\langle\xi^{*}(t)\nu(t^{\prime})\rangle$, and $\langle\nu^{*}(t)\nu(t^{\prime})\rangle$ may be modified without altering the expectation value of samples propagated by eq. (7). However, the generic constraint of non-negative probability does apply and sets limits to the range of choices. Consequently, the solutions of eq. (7) are of the form $\tilde{\rho}(t)=R_{L}(t)\rho_{S}(0)R^{\dagger}_{R}(t)$, where $R_{L/R}(t)$ is almost surely non-unitary.

As a consequence of non-unitarity, the variance of samples, expressed by observables (or a suitable metric applied to $\tilde{\rho}$), tends to increase with the length of the propagation time interval. The typical asymptotic behavior, shared with many processes involving multiplicative noise gardi09 , is exponential growth as a function of the interval length. Key quantities in this context are the noise spectra of $\nu^{\prime\prime}$ and $\xi^{\prime\prime}$, which drive both the sample trace and the Frobenius norm of $\tilde{\rho}$ towards a log-normal distribution. In some parameter regimes, this is harmless, in others, a finite-memory approach can halt this variance growth at a finite time stock16a . However, there is a problem in the strong-coupling limit, where the asymptotic dependence of sample variance on the coupling constant is exponential.

3 Optimization of stochastic driving

Varying the non-physical correlations with the aim of minimizing $\nu^{\prime\prime}$ and $\xi^{\prime\prime}$ seems a natural desideratum. Formally, this is an optimization problem in a function space with explicit equality constraints and inequality constraints implied by the non-negativity of probability measures.

With the notable exception of Ref. imai15 , this optimization was not attempted in previous work, and the non-physical correlations were chosen with an eye to easy numerical noise generation stock02 ; stock04 ; koch08 . $\xi$ was decomposed into a sum of independent terms $\xi^{l}\in\mathbbm{R}$ and $\xi^{s}\in\mathbbm{C}$ with

[TABLE]

where the prefactor $\mu$ of the Dirac delta function is chosen such that the time integral over (9) vanishes counterterm . All other correlations are zero, except $\langle\nu^{*}(t)\nu(t^{\prime})\rangle$ and $\langle(\xi^{(s)})^{*}(t)\xi^{(s)}(t^{\prime})\rangle$, which are fixed by setting them equal to each other and assigning them the minimal autocorrelation noise power allowed by eq. (9).

When minimizing the spectral noise power of $\nu^{\prime\prime}$ and $\xi^{\prime\prime}$, one needs to consider the full matrix of correlation functions after a Fourier transform,

[TABLE]

The rows and columns of $\mathsf{\Sigma}$ will be naturally labeled by indices $\xi^{\prime}$, $\xi^{\prime\prime}$, $\nu^{\prime}$ and $\nu^{\prime\prime}$. With the notations $S(\omega)=F[\langle\xi(t)\xi(t^{\prime})]$ and $D(\omega)=-iF[\langle\xi(t)\nu(t^{\prime})]$, the operator-noise identifications (4)–(6) now read

[TABLE]

Because the spectra of real-valued correlation functions are considered at this point, pairs $(-\omega,\omega)$ must be considered in the optimization. As a consequence, condition (5) leads to two separate equations (12) and (13).

The condition of $\mathsf{\Sigma}$ being positive semidefinite cannot easily be stated in the form of a simple algebraic expression. However, optimization problems involving semidefinite Hermitian matrices and linear constraints are special cases of convex optimization. These can be solved efficiently and reliably with modern numerical methods, implemented in packages such as CVX grant14 ; grant08 , which we used to minimize $\Sigma_{\xi^{\prime\prime}\xi^{\prime\prime}}+\Sigma_{\nu^{\prime\prime}\nu^{\prime\prime}}$, subject to all of the constraints mentioned above.

The numerical solutions invariably show specific features which guide us towards an analytic solution:

$\mathsf{\Sigma}$ is block diagonal; the pairs $(\xi^{\prime},\nu^{\prime\prime})$ and $(\xi^{\prime\prime},\nu^{\prime})$ are independent. 2. 2.

Each block is a rank 1 matrix.

With these provisos, the solution of the optimization problem can be written in closed form,

[TABLE]

with

[TABLE]

For an Ohmic reservoir, the resulting noise powers $\Sigma_{\xi^{\prime\prime}\xi^{\prime\prime}}$ and $\Sigma_{\nu^{\prime\prime}\nu^{\prime\prime}}$ are always smaller than in the previously used construction using (8) and (9) by at least 30%. The most significant advantage occurs in the case $|D|<S$, where the optimized values are $\Sigma_{\nu^{\prime\prime}\nu^{\prime\prime}}\approx|D|^{2}/S$ and $\Sigma_{\xi^{\prime\prime}\xi^{\prime\prime}}\approx|D|^{4}/S^{3}$, while the old approach had $\Sigma_{\xi^{\prime\prime}\xi^{\prime\prime}}\approx\Sigma_{\nu^{\prime\prime}\nu^{\prime\prime}}\approx|D|/2$. In the limit $|D|/S\to 0$, the optimal result coincides with the ansatz $\Sigma_{\xi^{\prime\prime}\xi^{\prime\prime}}=0$ made by Imai et al. for certain frequency ranges imai15 .

In the context of Ohmic dissipation at finite temperature, the low-frequency behavior is quite relevant: It is a most important parameter determining the growth of the sampling variance discussed above. At frequencies below the thermal energy, the noise powers of $\xi^{\prime\prime}$ and $\nu^{\prime\prime}$ are now smaller by a factor of $\beta\omega$, i.e., the noise powers $\Sigma_{\xi^{\prime\prime}\xi^{\prime\prime}}$ and $\Sigma_{\nu^{\prime\prime}\nu^{\prime\prime}}$ now vanish quadratically in the infrared (instead of linearly). The comparison shown in Figures 1 and 2 demonstrates the extremely beneficial effect of this on the convergence of simulation data, using the spin-boson system legge87 as an example. The symmetric spin-boson system is defined through $H_{S}=-\frac{\Delta}{2}\sigma_{x}$ and $A=\sigma_{z}$, using Pauli spin matrices. The reservoir is characterized by an Ohmic spectral density parameterized by UV cutoff frequency $\omega_{c}$ and a dimensionless dissipation constant $K$ weiss08b ; stock04 . Figure 1 shows simulation results with identical parameters, including number of samples using the optimimzed noise statistics, eq. (23), represented by smoother lines and the conventional statistics (8)–(9), represented by lines indentifyable by visible growth of their noise amplitude.

Fig. 2 and table 2 compare the performance characteristics of the old and new methods in further detail. Depending on temperature, the new approach reduces the growth rate of the sample variance significantly. The resulting savings factor in the required number of samples is therefore also an exponential function of simulation time. Savings factors of up to three orders of magnitude allow simulations extending up to approximate equilibration, using optimized samples. The savings factor is largest for high temperature and strong coupling. However, the asymptotic numerical cost still grows exponentially in the limit of very strong coupling.

4 Conclusions

The noise samples used in an SLN-based propagation of open quantum systems can be optimized towards small imaginary parts. The resulting reduction in the number of required samples can be several orders of magnitude. No additional numerical costs arise from the optimization since it has an analytic solution. The range of applications for SLN simulations thus widens in virtually all parameter regimes, moderate or high temperature being the regime with the greatest benefit. With increasing availability of parallel computing resources, SLN-based methods are becoming an attractive, versatile tool for the dynamics of open quantum systems.

The authors would like to thank two anonymous referees for valuable suggestions. K.S. showed that standard methods of semidefinite programming apply to the noise optimization and obtained the numerical optimization results; J.T.S. derived the analytic result (23). Both authors contributed equally to the rest of the material.