Optimized auxiliary oscillators for the simulation of general open quantum systems

TL;DR

This paper introduces a systematic method to construct few-body damped harmonic oscillator networks that accurately simulate general bosonic environments in open quantum systems, enabling efficient and flexible modeling of non-Markovian effects.

Contribution

The paper presents a novel algorithm for designing effective oscillator networks that replicate complex environments with controlled error, applicable to a wide range of open quantum systems.

Findings

Accurately reproduces environment effects with few oscillators

Enables efficient simulation of non-Markovian dynamics

Applicable at arbitrary temperatures, including zero

Abstract

A method for the systematic construction of few-body damped harmonic oscillator networks accurately reproducing the effect of general bosonic environments in open quantum systems is presented. Under the sole assumptions of a Gaussian environment and regardless of the system coupled to it, an algorithm to determine the parameters of an equivalent set of interacting damped oscillators obeying a Markovian quantum master equation is introduced. By choosing a suitable coupling to the system and minimizing an appropriate distance between the two-time correlation function of this effective bath and that of the target environment, the error induced in the reduced dynamics of the system is brought under rigorous control. The interactions among the effective modes provide remarkable flexibility in replicating non-Markovian effects on the system even with a small number of oscillators, and the resulting Lindblad equation may therefore be integrated at a very reasonable computational cost using standard methods for Markovian problems, even in strongly non-perturbative coupling regimes and at arbitrary temperatures including zero. We apply the method to an exactly solvable problem in order to demonstrate its accuracy, and present a study based on current research in the context of coherent transport in biological aggregates as a more realistic example of its use; performance and versatility are highlighted, and theoretical and numerical advantages over existing methods, as well as possible future improvements, are discussed.

Tables

Table 1. Table 1: Effective parameters for N = 4 𝑁 4 N=4 surrogate modes corresponding to the correlation function of an Ohmic bath at temperature T = Ω c 𝑇 subscript Ω 𝑐 T=\Omega_{c} (Eq. ( 27 )). All parameters have dimensions of frequency and are given in units of Ω c subscript Ω 𝑐 \Omega_{c} ; the last c n subscript 𝑐 𝑛 c_{n} is real.

	Mode 1		Mode 2		Mode 3		Mode 4
${\Omega }_n$	$0.512683$		$2.53779$		$4.53293$		$0.151433$
$g_n$		$1.82454$		$3.20774$		$1.60194$
${\Gamma }_n$	$0.056336$		$4.42709$		$15.7371$		$0.110104$
$c_n$	$-0.962917$		$-0.227707$		$0.231179$		$0.818093$
	$+0.819128i$		$+0.0701249i$		$-0.137866i$

Table 2. Table 2: Total simulation times, absolute values of the dipole correlation function C μ ​ ( t ) subscript 𝐶 𝜇 𝑡 C_{\mu}(t) of the dimer at the final time and DAMPF bond dimensions (if applicable) for J 0 ​ ( ω ) subscript 𝐽 0 𝜔 J_{0}(\omega) and J 1 ​ ( ω ) subscript 𝐽 1 𝜔 J_{1}(\omega) , respectively. The simulations at T = 300 ​ K 𝑇 300 K T=300\,\mathrm{K} were performed with a time-adaptive bond dimension.

$T$ $\left(K\right)$	$t_{\max }$ $\left(\mathrm{ps}\right)$		$\left|C_{\mu }\left(t_{\max }\right)\right|$		$\chi$
	$J_0\left(\omega \right)$	$J_1\left(\omega \right)$	$J_0\left(\omega \right)$	$J_1\left(\omega \right)$
$0$	$20.0$	$11.0$			$12$
$77$	$6.75$	$5.5$		$\sim {10}^{-3}$	$12$
$300$	$0.27$	$0.17$	$\sim {10}^{-3}$	$\sim {10}^{-3}$	variable

Table 3. Table 4: Ohmic spectral density with cutoff frequency Ω c subscript Ω 𝑐 \Omega_{c} , temperature T = Ω c 𝑇 subscript Ω 𝑐 T=\Omega_{c} : parameters in units Ω c subscript Ω 𝑐 \Omega_{c} and local dimensions.

	Mode 1		Mode 2		Mode 3		Mode 4
${\Omega }_n$	$0.512683$		$2.53779$		$4.53293$		$0.151433$
$g_n$		$1.82454$		$3.20774$		$1.60194$
${\Gamma }_n$	$0.056336$		$4.42709$		$15.7371$		$0.110104$
$c_n$	$-0.962917$		$-0.227707$		$0.231179$		$0.818093$
	$+0.819128i$		$+0.0701249i$		$-0.137866i$
$d_{\mathrm{loc}}$ (spin)	$5$		$4$		$4$		$7$
$d_{\mathrm{loc}}$ (chain)	$7$		$4$		$3$		$8$

Table 4. Table 5: Ohmic spectral density with cutoff frequency Ω c subscript Ω 𝑐 \Omega_{c} , temperature T = 5 2 ​ Ω c 𝑇 5 2 subscript Ω 𝑐 T=\frac{5}{2}\Omega_{c} : parameters in units Ω c subscript Ω 𝑐 \Omega_{c} and local dimensions.

	Mode 1		Mode 2		Mode 3		Mode 4		Mode 5
${\Omega }_n$	$0.306859$		$0.361308$		$0.167597$		$0.0297981$		$0.00236395$
$g_n$		$4.17718$		$2.1243$		$0.673391$		$0.166947$
${\Gamma }_n$	$16.0093$		$2.76375$		$0.00358704$		$0.0949691$		$0.0517414$
$c_n$	$-0.166675$		$0.21927$		$1.61933$		$0.187388$		$1.1553$
	$-0.0342019i$		$+0.103791i$		$-0.00703994i$		$-1.07416i$
$d_{\mathrm{loc}}$	$3$		$3$		$4$		$4$		$6$

Table 5. Table 7: Adolphs–Renger spectral density, temperature T = 77 ​ K 𝑇 77 K T=77\,\mathrm{K} : parameters in units u = 100 ​ cm − 1 𝑢 100 superscript cm 1 u=100\,\mathrm{cm}^{-1} and local dimensions.

	Mode 1		Mode 2		Mode 3		Mode 4
${\Omega }_n$	$3.05106$		$2.74196$		$0.00670418$		$0.00780109$
$g_n$		$2.74161$		$2.01796$		$0.33975$
${\Gamma }_n$	$0.0284151$		$11.6481$		$0.00549033$		$0.0184315$
$c_n$	$-0.910465$		$-0.135049$		$0.524001$		$0.114767$
	$-0.0164266i$		$-0.0104797i$		$+0.317767i$
$d_{\mathrm{loc}}$	$5$		$4$		$6$		$8$

Table 6. Table 8: Adolphs–Renger spectral density, temperature T = 300 ​ K 𝑇 300 K T=300\,\mathrm{K} : parameters in units u = 100 ​ cm − 1 𝑢 100 superscript cm 1 u=100\,\mathrm{cm}^{-1} and local dimensions.

	Mode 1		Mode 2		Mode 3		Mode 4
${\Omega }_n$	$0.788783$		$0.414407$		$-0.0300357$		$-0.034035$
$g_n$		$3.10576$		$0.978945$		$0.294823$
${\Gamma }_n$	$10.4575$		$0.0934767$		$0.00983292$		$0.0167273$
$c_n$	$0.189405$		$1.23326$		$0.0221509$		$0.365249$
	$+0.0639657i$		$+0.451035i$		$+0.962709i$
$d_{\mathrm{loc}}$	$3$		$4$		$7$		$7$

Table 7. Table 10: Antisymmetrized Lorentzian spectral density with Ω = 200 ​ cm − 1 Ω 200 superscript cm 1 \Omega=200\,\mathrm{cm}^{-1} , Γ = 10 ​ cm − 1 Γ 10 superscript cm 1 \Gamma=10\,\mathrm{cm}^{-1} , S = 0.25 𝑆 0.25 S=0.25 : parameters in units u = 100 ​ cm − 1 𝑢 100 superscript cm 1 u=100\,\mathrm{cm}^{-1} and local dimensions.

	$T=0K$	$T=77K$			$T=300K$
		Mode 1		Mode 2	Mode 1		Mode 2
${\Omega }_n$	$2.00$	$-0.318699$		$0.316331$	$-0.00048954$		$0.000480821$
$g_n$	$-$		$1.976$			$2.00052$
${\Gamma }_n$	$0.098296$	$0.045988$		$0.138442$	$0.00953908$		$0.190362$
$c_n$	$0.992322$	$0.764199$		$0.676024$	$1.45733$		$0.343374$
		$+0.000002i$			$+0.000003i$
$d_{\mathrm{loc}}$	$6$	$5$		$6$	$8$		$8$

Equations

$$H:=H_S\otimes {\mathbb{I}}_E+{\mathbb{I}}_S\otimes H_E+H_I$$

$$H_E:={\int }_0^{\infty }\mathrm{d}\omega \hslash \omega a_{\omega }^{†}{a}_{\omega }$$

$$H_I:=\sum _k{A}_{Sk}\otimes G_{Ek}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)=-\frac{i}{\hslash }[H,\rho (t)]$$

$${\rho }_S(t)={\mathrm{Tr}}_E[\rho (t)].$$

$$⟨G_{Ek}(t){⟩}_E$$

$$C_{k{k}^{\prime }}^E(t+\tau ,\tau )$$

$$H^{\prime }:=H_S\otimes {\mathbb{I}}_R+{\mathbb{I}}_S\otimes H_R+H_I^{\prime }$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_R(t)={\mathcal{L}}_R[{\rho }_R(t)],$$

$${\mathcal{L}}_R[{\rho }_R]:=-\frac{i}{\hslash }[H_R,{\rho }_R]+{\mathcal{D}}_R[{\rho }_R]$$

$${\mathcal{D}}_R[{\rho }_R]:=\sum _{i,j=1}^m{\Gamma }_{ij}\left(L_{Ri}{\rho }_R{L}_{Rj}^{†}-\frac{1}{2}\left\{L_{Rj}^{†}{L}_{Ri},{\rho }_R\right\}\right)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)=\mathcal{L}[\rho (t)]$$

$$\mathcal{L}[\rho ]:=-\frac{i}{\hslash }[H^{\prime },\rho ]+\mathcal{D}[\rho ]$$

$$H_R:=\sum _n\hslash {\omega }_n{b}_n^{†}{b}_n$$

$$H_I^{\prime }:=\sum _l{A}_{Sl}\otimes F_{Rl}$$

$$⟨F_{Rl}(t){⟩}_R$$

$$C_{l{l}^{\prime }}^R(t+\tau ,\tau )$$

$${\rho }_S^U(t):={\mathrm{Tr}}_E[\rho (t)]$$

$${\rho }_S^L(t):={\mathrm{Tr}}_R[\rho (t)]$$

$$⟨F_{Rk}(t){⟩}_R=⟨G_{Ek}(t){⟩}_E\forall k,t$$

$$C_{k{k}^{\prime }}^R(t+\tau ,\tau )=C_{k{k}^{\prime }}^E(t+\tau ,\tau )\forall k,k^{\prime },t,\tau ,$$

$${\rho }_S^L(t)={\rho }_S^U(t)\forall t.$$

$$C_{k{k}^{\prime }}^E(t+\tau ,\tau )=C_{k{k}^{\prime }}^E(t,0),$$

$$F_{Rk}=\sum _n(c_{nk}{b}_n+c_{nk}^{\ast }{b}_n^{†}).$$

$${\mathcal{L}}_R[{\rho }_{0R}]=0.$$

$$H_R:=\sum _{n=1}^N\hslash {\Omega }_n{b}_n^{†}{b}_n+\sum _{n=1}^{N-1}\hslash g_n\left(b_n{b}_{n+1}^{†}+b_n^{†}{b}_{n+1}\right),$$

$$\begin{gathered} \frac{\mathrm{d}}{\mathrm{d}t}{\rho }_R(t)=-\frac{i}{\hslash }[H_R,{\rho }_R(t)] \\ +\sum _{n=1}^N{\Gamma }_n\left(b_n{\rho }_R(t)b_n^{†}-\frac{1}{2}\left\{b_n^{†}{b}_n,{\rho }_R(t)\right\}\right) \end{gathered}$$

$${\rho }_{0R}=\underset{n=1}{\overset{N}{⨂}}|0⟩{⟨0|}_n.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}⟨b_n(t){⟩}_R=\sum _{m=1}^N{M}_{nm}⟨b_m(t){⟩}_R,$$

$$M_{nm}:=\left(-\frac{{\Gamma }_n}{2}-i{\Omega }_n\right){\delta }_{nm}-i(g_m{\delta }_{nm+1}+g_{m-1}{\delta }_{nm-1}).$$

$$C_{k{k}^{\prime }}^R(t)=\sum _{n=1}^N(w_n{)}_{k{k}^{\prime }}{e}^{{\lambda }_{n}t}$$

Full text

I Introduction

Any physical system in nature may be studied theoretically in complete isolation from its surroundings. However, since interactions with uncontrolled environmental degrees of freedom are unavoidable in practice, this condition is never actually realized. The effects of said degrees of freedom on the dynamics and general properties of a system are especially important in quantum mechanics, where the time and energy scales involved are likely to make interactions between the system and the surrounding environment a key actor in their own right in the physics at play. The goal of the theory of open quantum systems is to determine the behavior and investigate the physical properties of systems both in and out of equilibrium by properly accounting for environmental effects and other external influences (e.g. driving forces) using appropriate analytical or numerical methods (BreuerPetruccione, ; GardinerZoller, ; Weiss, ; RivasHuelga, ).

The starting point of such methods may be either a microscopic model for the system and the environment, such as the spin-boson (Leggett_SpinBoson, ; BreuerPetruccione, ), Caldeira–Leggett (CaldeiraLeggettModel, ) or more complex models, or an effective description of the system alone with the effects of the environment implicitly taken into account via a quantum master equation (GoriniKossakowskiSudarshan, ; GoriniFrigerioVerriKossakowskiSudarshan, ; Lindblad, ; Nakajima, ; Zwanzig, ; Prigogine, ; Shibata_TCL, ). The former setup leads to a wide variety of potentially more complete and general treatments, but this greater range of attainable results and predictions comes at moderate to high computational costs (NEGF_Review, ; Tanimura_HEOM, ; Makri_QUAPIletter, ; ChinPlenio_TEDOPA, ; PriorPlenio_TEDOPA, ; DiosiStrunz_QuantumStateDiffusion, ; Piilo_NonMarkovianQuantumJumpsPRL, ); the latter construction is typically less expensive but applies to a constrained class of physical settings, since it either delivers accurate results only in a few well-defined limiting cases (GoriniFrigerioVerriKossakowskiSudarshan, ; Davies_WeakCoupling, ; DumckeSpohn_WeakCoupling, ) or relies on equations which are difficult to derive for general systems (Nakajima, ; Zwanzig, ; Prigogine, ; Shibata_TCL, ; SmirneVacchini_NakajimaVsTCL, ; BreuerKapplerPetruccione_TCL, ). Provided that the necessary assumptions on the system-environment interaction are satisfied, however, efficient methods for the solution of the master equation are widely available (Minchev_MatExp, ; PlenioKnight_MCWF, ; GisinPercival_QuantumStateDiffusion, ).

Much theoretical research in recent decades has focused on the study of complex non-Markovian environments (RivasHuelgaPlenio_NonMarkovianity, ; BreuerVacchini_NonMarkovianity, ; DeVegaAlonso_NonMarkovianity, ; LiHallWiseman_NonMarkovianity, ), for which analytical results are hard to obtain except for specific models, and numerical simulation may become very challenging depending on the physical regime of interest. For thermal bosonic environments, the most commonly studied category by far, numerical methods developed for a general treatment of non-Markovian problems include e.g. Hierarchical Equations of Motion (HEOM) (Tanimura_HEOM, ; TanimuraKubo_HEOM, ), Quasi-Adiabatic Path Integrals (QUAPI) (Makri_QUAPIletter, ; TopalerMakri_QUAPI, ; MakriMakarov_QUAPI1, ; MakriMakarov_QUAPI2, ), Nonequilibrium Green’s Function (NEGF) techniques (Danielewicz_NEGF, ; NEGF_Review, ), Non-Markovian Quantum State Diffusion (NMQSD) and similar stochastic methods (Diosi_QuantumStateDiffusion, ; Strunz_QuantumStateDiffusion, ; DiosiStrunz_QuantumStateDiffusion, ; Piilo_NonMarkovianQuantumJumpsPRA, ; Piilo_NonMarkovianQuantumJumpsPRL, ), Time-Evolving Matrix Product Operators (TEMPO) (Lovett_TEMPO, ) or simulated evolution of the state using the time-adaptive Density Matrix Renormalization Group (t-DMRG) (Daley_tDMRG, ; Guifre_tDMRG, ; Schollwock_tDMRG, ) in combination with convenient exact mappings of the environment e.g. into one-dimensional oscillator chains well suited for such numerical methods, as in the Time-Evolving Density with Orthogonal Polynomials Algorithm (TEDOPA) (ChinPlenio_TEDOPA, ; PriorPlenio_TEDOPA, ; TamascelliSmirne_ThermalizedTEDOPA, ), to name a few. These methods are often referred to as numerically exact, in the sense that they are designed to address problems from the bottom up, requiring only numerical approximations (e.g. Hilbert space truncation, discretized integrals or finite expansions of relevant functions) in order to keep the costs manageable, but otherwise posing no physical restrictions on the models themselves; these numerical errors can sometimes be bounded rigorously, e.g. for TEDOPA (WoodsCramerPleino_TEDOPAerrorbars, ; WoodsPlenio_LiebRobinsonBounds, ) or HEOM (SpinBosonBounds, ). Finite bosonic environments (PiiloManiscalco_HarmonicNetworks, ) can also be used as an approximate treatment for simulation times short enough to prevent recurrence in the dynamics.

An alternative route for the numerical study of such nontrivial open-system problems is to model environmental effects on a system by splitting them into coherent, information-preserving contributions and purely dissipative Markovian damping. Then one can devise effective models in which the system of interest is coupled explicitly to a finite auxiliary system acting as the non-Markovian core of the environment, and dissipation is accounted for through Markovian damping of these auxiliary degrees of freedom. This is the idea underlying approaches such as the pseudomode method (Imamoglu_pseudomodes, ; Garraway_pseudomodes, ; Dalton_pseudomodes, ; Lemmer_SpinBoson, ), the reaction-coordinate mapping (IlesSmith_ReactionCoordinate, ; IlesSmith_StructuredEnvironments, ; Lambert_SpinBosonStudy, ) or other techniques based on the same concept but differing in the ansatz used to create the effective environment and the techniques to solve for the dynamics (Fruchtman_Perturbative, ; Falci_1overf, ; Schwarz_LindbladDrivenDiscretizedLeads, ; Luchnikov_TensorNetworkLindblad, ; Faccioli_QTFT, ). Such remappings of open-system problems can be very convenient numerically, but are not always grounded in a mathematically rigorous and physically sound relation between the original and effective environments, making assessment of their accuracy somewhat challenging.

In this paper, we present a new approach to general open quantum systems interacting with Gaussian bosonic environments. Our method combines the simplicity and efficiency of simulating a small set of effective degrees of freedom with analytical equivalence relations between the structure and parameters of this auxiliary system and the exact properties of the microscopic environment. Even in cases in which no exact equivalence holds, the physical error from replacing a unitary environment by a dissipative one is kept to a bare minimum and under rigorous control.

Our scheme is based on a quantitatively certified recipe to construct networks of interacting, damped harmonic oscillators specifically designed to mimic any given target environment as specified by its spectral density and temperature. The reduced dynamics is then computed by solving a time-homogeneous quantum master equation of the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) type (GoriniKossakowskiSudarshan, ; GoriniFrigerioVerriKossakowskiSudarshan, ; Lindblad, ) for the system coupled to these effective harmonic modes and tracing them out at the end. The theoretical foundation underlying this construction lies in a recently proved equivalence theorem between unitary and non-unitary Gaussian environments in open quantum systems (TSO_Theorem, ), which states that the reduced dynamics of a system coupled to an environment of either type is identical if the single-time averages and two-time correlation functions of the environment operators relevant to the interaction are the same. Exploiting this notion, we introduce a systematic procedure by which the effective environment is tailored to reproduce the correlation function of the target environment with an accuracy quantitatively controlled through known error bounds for Gaussian environments (SpinBosonBounds, ). The advantages of the method proposed are the simple yet versatile structure of the effective environments, which can emulate a broad range of nontrivial unitary environments using small numbers of auxiliary modes, the small, controlled error in the resulting effective dynamics, a high flexibility in the physical regimes which can be studied at comparably low costs, such as high and low temperature and strong as well as weak coupling, and numerical simplicity, since the simulations only require solving a Lindblad equation.

We have organized the presentation of our results as follows: in Section II we will outline the theoretical background and state the equivalence theorem from Ref. (TSO_Theorem, ) lying at the core of our method; Section III details the procedure by which an effective environment corresponding to a nontrivial microscopic one may be constructed, and includes an analysis of the theoretical implications and approximations involved; a demonstration of our scheme on the spin-boson model as an exactly solvable test system, with accuracy and performance reports as well as a profile of the numerical advantages and disadvantages of the method in different physical regimes, is given in Section IV; Section V contains two applications of the method to systems in structured environments relevant to current research, namely optical signatures of coherent effects in biomolecular aggregates and the propagation of excitations in organic polymers with photovoltaic properties; in Section VI we discuss the current state of the method, focusing on its scope and applicability, numerical and conceptual strengths and limitations as well as some possible improvements; finally, Section VII summarizes our conclusions and future prospects.

II Theoretical foundations and scope of the method

The non-perturbative method we are going to introduce relies on the equivalence theorem between unitary and dissipative environments stated and proved in Ref. (TSO_Theorem, ); in order to set the stage for discussing our work, we will now introduce the relevant notation, outline the physical context in which the theorem applies and state it explicitly for reference within the paper.

II.1 Gaussian unitary environments

A wide array of open quantum system (OQS) problems, ranging e.g. from quantum Brownian motion (BreuerPetruccione, ; CaldeiraLeggettModel, ) to dissipative cavity and circuit electrodynamics (YurkeDenker_QuantumCircuits, ; RibeiroVieira_Transport, ) or the study of charge and energy transfer in noisy natural or artificial aggregates (Scholes_QBioNature, ; HuelgaPlenio_QBio, ), can be modeled microscopically by coupling the system of interest to an infinite collection of harmonic oscillators: the full system-environment Hamiltonian takes the form

[TABLE]

where $H_{S}$ is the free Hamiltonian of the system,

[TABLE]

is the free Hamiltonian of the environment, expressed in terms of creation and annihilation operators obeying the continuum canonical commutation relations $[a_{\omega},a_{\omega^{\prime}}]=0$, $[a_{\omega},a^{\dagger}_{\omega^{\prime}}]=\delta(\omega-\omega^{\prime})$, and the two are coupled through a general interaction term of the form (BreuerPetruccione, )

[TABLE]

with $A_{Sk}$ and $G_{Ek}$ operators acting on the system and the environment, respectively. In the following, we will consider these to be Hermitian, without loss of generality (RivasHuelga, ).

The global state $\rho$ of the system and the environment at time $t$ is determined by the Liouville–Von Neumann equation

[TABLE]

and the initial state $\rho_{0}\coloneqq\rho(0)$; the reduced state $\rho_{S}$ of the system at time $t$ is obtained by taking the partial trace over the environmental degrees of freedom:

[TABLE]

We are interested in the reduced dynamics of systems interacting with Gaussian environments, i.e. with $H_{I}$ linear in $a_{\omega}$ and $a^{\dagger}_{\omega}$ and factorizing initial conditions $\rho_{0}=\rho_{0S}\otimes\rho_{0E}$ with $\rho_{0E}$ a Gaussian state. Then the oscillators can be traced out exactly, and the reduced dynamics of the system only depends on the single- and two-time environmental averages $\langle G_{Ek}(t)\rangle_{E}$ and $C^{E}_{kk^{\prime}}(t+\tau,\tau)\coloneqq\langle G_{Ek}(t+\tau)G_{Ek^{\prime}}(\tau)\rangle_{E}$ as given by the evolution of the oscillators with no coupling to the system:

[TABLE]

with $U_{E}(t)\coloneqq e^{-iH_{E}t/\hbar}$.

II.2 Gaussian dissipative environments

Considering infinite environments evolving unitarily in the absence of a coupled system is one way to bring about dissipation and decoherence in the evolution of the latter when the coupling is nonzero. Alternatively, one may consider finite environments which evolve non-unitarily according to a quantum master equation (QME). In this case, one may start from a different combined Hamiltonian

[TABLE]

and a QME describing the evolution of the dissipative environment when decoupled from the system:

[TABLE]

where the new quantum Liouville superoperator

[TABLE]

on the right-hand side includes a dissipator

[TABLE]

with a positive semidefinite rate matrix $\Gamma_{ij}$. This makes the dynamics non-unitary but ensures a completely positive and trace-preserving evolution at all positive times; the rate matrix $\Gamma_{ij}$, which we take to be constant, can always be brought into diagonal form by changing the basis of operators $L_{Ri}$ (BreuerPetruccione, ), giving the quantum dynamical semigroup master equation for Markovian open systems first derived by Gorini, Kossakowski, Sudarshan and Lindblad (GoriniKossakowskiSudarshan, ; GoriniFrigerioVerriKossakowskiSudarshan, ; Lindblad, ). We will refer to this master equation simply as the Lindblad equation throughout this paper.

The full state of the system and a non-unitary environment evolves according to the QME

[TABLE]

where

[TABLE]

is the complete quantum Liouvillian for the system and the environment and $\mathcal{D}\coloneqq\mathbb{I}\otimes\mathcal{D}_{R}$ embeds the dissipator $\mathcal{D}_{R}$ into the full Hilbert space of the problem.

For harmonic environments coupled linearly to the system, i.e. for

[TABLE]

with $[b_{m},b_{n}]=0$, $[b_{m},b^{\dagger}_{n}]=\delta_{mn}$ and

[TABLE]

with $F_{Rl}$ linear in the creation and annihilation operators, if one also takes the Lindblad operators $L_{Ri}$ linear in $b_{n}$ and $b^{\dagger}_{n}$ and initial conditions $\rho_{0}=\rho_{0S}\otimes\rho_{0R}$ with a Gaussian $\rho_{0R}$, then the reduced dynamics of the system will only depend on the environment through $\langle F_{Rl}(t)\rangle_{R}$ and $C^{R}_{ll^{\prime}}(t+\tau,\tau)\coloneqq\langle F_{Rl}(t+\tau)F_{Rl^{\prime}}(\tau)\rangle_{R}$, again considering the free dynamics of the environment with no system attached, like in the unitary case:

[TABLE]

Note that the two-time correlation function (10) has the form one would obtain by applying the quantum regression hypothesis (Lax_QRT, ), which must be handled with some care in general but is true by construction for the Lindblad-damped environments relevant to our work. No approximation is required or implied at this stage (TSO_Theorem, ).

II.3 Equivalence between unitary and non-unitary environments

While it is clear that if two unitary Gaussian environments share the same averages $\langle G_{Ek}(t)\rangle_{E}$ and correlation functions $C^{E}_{kk^{\prime}}(t+\tau,\tau)$ at all times they will give rise to the same reduced dynamics if coupled to a system, this is not obvious if one or both environments are not unitary. In Ref. (TSO_Theorem, ) it was shown, using the unitary dilation formalism for Lindblad equations (GardinerZoller, ), that this still holds for non-unitary environments under the same conditions. We restate this result here for reference.

Define the reduced dynamics

[TABLE]

for some system $S$ coupled to a unitary environment and evolving according to Eq. (2) from factorizing initial conditions with the environment starting in a Gaussian state, and the reduced dynamics

[TABLE]

for the same system coupled to a non-unitary environment and evolving according to Eq. (8) from factorizing initial conditions with the environment starting in a Gaussian state.

Both environments are taken to be harmonic and coupled to the system through the same set of $A_{Sk}$ operators in $H_{I}$ and $H^{\prime}_{I}$, with the corresponding $G_{Ek}$ and $F_{Rk}$ as well as the Lindblad operators $L_{Ri}$ of the non-unitary environment linear in the relevant creation and annihilation operators. The initial state of the system is taken to be the same.

Theorem 1

(TSO_Theorem*, *)** Under the above assumptions, if

[TABLE]

and

[TABLE]

then

[TABLE]

This theorem is the cornerstone of our method; for the sake of clarity and an easier understanding of the rest of this paper, some remarks are in order.

First of all, it is important to stress that Gaussianity is a key ingredient of Theorem 1, because in principle all correlation functions up to infinite order would have to be equal for two environments to have the same effect on a system, but for Gaussian environments the single- and two-time functions generate all the others. This restricts the initial state of the environment to the Gaussian family; in this work, we will only consider system-environment product states which are Gaussian in the environmental variables as initial states, leaving the free Hamiltonian, interaction operators and initial density matrix of the system arbitrary.

Furthermore, we restrict our study to systems coupled to bosonic baths in this paper but a result analogous to Theorem 1 was recently proved for fermionic environments as well (Chen_FermionicTSO, ), making an extension of our work to fermionic open-system problems possible.

Finally, for physical reasons discussed in Section VI and thoroughly analyzed in Ref. (Talkner_NoQRT, ), in general a finite network of damped harmonic oscillators does not yield a two-time correlation function exactly equal to that of an infinite bath in a Gaussian equilibrium state, so we will apply the theorem in approximate form by looking for effective parameters such that $C^{R}_{kk^{\prime}}(t+\tau,\tau)\approx C^{E}_{kk^{\prime}}(t+\tau,\tau)$ and hence $\rho^{L}_{S}(t)\approx\rho^{U}_{S}(t)$ (single-time expectation values of coupling operators are typically zero or can be set to zero and will no longer be dealt with in this work), relying on the fact that the error in the former approximate relation rigorously bounds that in the latter, as established in previous work (SpinBosonBounds, ).

Other than these caveats, no further problems arise in terms of applicability; in particular, temperature and coupling strength between system and environment pose no theoretical or computational limits in principle.

In the next sections we will show how one may exploit the theorem to systematically construct simple networks of damped harmonic oscillators, which can stand in for complex, highly non-Markovian thermal baths at any temperature, by comparing the associated correlation functions (10) and (5). This procedure is independent of the system and the effective environments obtained through it can then be coupled arbitrarily strongly to any system of interest. Standard simulation methods for Lindblad equations may then be used to obtain the reduced dynamics at potentially very low computational costs.

III Systematic construction of effective environments

From now on, we will consider unitary environments with Gaussian stationary states, such as thermal baths, and assume them to be initialized in such states, so that

[TABLE]

which will be denoted by $C^{E}_{kk^{\prime}}(t)$ in the following.

Any harmonic oscillator network obeying a Lindblad equation of the form (7), with $\mathcal{L}_{R}$ quadratic in $b_{n}$ and $b^{\dagger}_{n}$, must also start from a stationary $\rho_{0R}$ in order to give a time-homogeneous correlation function matrix $C^{R}_{kk^{\prime}}(t)\coloneqq C^{R}_{kk^{\prime}}(t+\tau,\tau)=C^{R}_{kk^{\prime}}(t,0)$ for operators of the form

[TABLE]

The condition for $\rho_{0R}$ to be a stationary state of Eq. (7) is

[TABLE]

For the initial state of our effective environments, we will therefore need a Gaussian $\rho_{0R}$ satisfying this property.

III.1 Ansatz and correlation function structure

The correlation functions $C^{R}_{kk^{\prime}}(t)$ of the auxiliary environment depend on all parameters appearing in $\mathcal{L}_{R}$, $\rho_{0R}$ and the operators $F_{Rk}$: unrestricted geometries and initial states allow for more generality at the expense of keeping potentially redundant parameters in the model and restricting the range of properties that can be easily calculated; to strike a balance between simplicity and versatility, we will now take an ansatz for the configuration and initial density matrix of the surrogate oscillator network such that the quantities of interest have a simple expression with little loss of generality; for a more extensive discussion of the technical details, we refer the reader to Appendix A.

We choose a free Hamiltonian $H_{R}$ corresponding to a chain of $N$ oscillators with a hopping interaction between nearest neighbors:

[TABLE]

where the couplings $g_{n}$, as well as one of the coefficients $c_{nk}$ in the interaction operators $F_{Rk}$ appearing in $C^{R}_{kk^{\prime}}(t)$, can be assumed real without loss of generality if the $F_{Rk}$ are nonlocal, i.e. acting on all effective modes (see Appendix A). We complete the QME by adding local thermal dissipators at zero temperature acting on each oscillator:

[TABLE]

so that the stationary initial state satisfying Eq. (14) is just the overall vacuum state

[TABLE]

Note that a zero-temperature master equation for the effective environment does not restrict the temperature of the target environments it can simulate; the effect of a nonzero temperature in the target bath will simply be encoded in the parameters of the oscillator network, as is done in other approaches (MayKuehn, ; DeVegaBanuls_Thermofield, ; DiosiGisinStrunz_NonMarkovianity, ; Ritschel_AbsorptionSpectra, ; TamascelliSmirne_ThermalizedTEDOPA, ).

The QME (16) and initial condition (17) lead to two decoupled sets of linear equations for $\langle b_{n}(t)\rangle_{R}$ and $\langle b^{\dagger}_{n}(t)\rangle_{R}$ related by Hermitian conjugation. For $\langle b_{n}(t)\rangle_{R}$ one has

[TABLE]

where

[TABLE]

The two-time correlation function $\langle b_{n}(t)b^{\dagger}_{m}(0)\rangle_{R}$ also evolves according to Eq. (18) as a direct consequence of the quantum regression hypothesis, which holds by construction in this context and states that correlation functions $\langle A(t+\tau)B(t)\rangle$ obey the same equations of motion as the single-time expectation values $\langle A(t+\tau)\rangle$ (GardinerZoller, ; Carmichael, ). This is equivalent to the statement that they can be written in the explicit form given in Eq. (10). Integrating Eq. (18) for $\langle b_{n}(t)b^{\dagger}_{m}(0)\rangle_{R}$ and plugging the result (as well as its conjugate $\langle b^{\dagger}_{n}(t)b_{m}(0)\rangle_{R}$, which is identically zero for our initial state (17)) into the expression of $C^{R}_{kk^{\prime}}(t)$ in terms of the operators $F_{Rk}$ as given in Eq. (13), one finds

[TABLE]

where $\lambda_{n}$ are the eigenvalues of the matrix $M$ defined in Eq. (19), which we assume to be non-degenerate for simplicity (see Appendix A for further discussions), and

[TABLE]

are complex coefficients obtained from the definition of the operators $F_{Rk}$ and the left and right eigenvectors $\mathbf{u}^{n}$ and $\mathbf{v}^{n}$ corresponding to each $\lambda_{n}$, normalized in such a way that $\sum^{N}_{l=1}v^{m}_{l}u^{n}_{l}=\delta_{mn}$. This exponential structure is a consequence of the Lindblad dynamics of the effective environment, which is a requirement of Theorem 1.

III.2 Transformation to Surrogate Oscillators

Consider an OQS problem described by a microscopic model of the form (1); for simplicity we will now assume a single interaction term, to which there corresponds a correlation function $C^{E}(t)$. Our goal is to find the matrix elements $M_{mn}$ and operator coefficients $c_{n}$ of some operator $F_{R}$ as given in Eq. (13) such that the resulting effective correlation function $C^{R}(t)$ is as close as possible to $C^{E}(t)$.

The form of $C^{R}(t)$ in terms of $M_{mn}$ and $c_{n}$ is given by Eq. (20), where the eigenvalues $\lambda_{n}$ and weights $w_{n}$ can be thought of as functions of the free parameters $\Omega_{k}$, $g_{k}$, $\Gamma_{k}$ and—only the $w_{n}$—$c_{k}$ with $n,k=1,\dots,N$, where $N$ is the number of oscillators making up the effective bath.

In order to determine the values of these free parameters such that $C^{R}(t)\approx C^{E}(t)$, we proceed in two steps. First, we perform a nonlinear fit on $C^{E}(t)$ using $N$ damped exponentials with complex coefficients

[TABLE]

for instance using Prony analysis (Marple_Prony, ). Note that the sum of the coefficients $\tilde{w}_{n}$ is real and positive, since $\sum^{N}_{n=1}\tilde{w}_{n}=\tilde{C}^{E}(0)$ must be equal to $C^{E}(0)=\left\langle G^{2}_{E}(0)\right\rangle_{E}$, which is positive because $G_{E}$ is Hermitian.

Then we solve the problem of matching or getting as close as possible to the $\tilde{\lambda}_{n}$ and $\tilde{w}_{n}$ with the $\lambda_{n}(\Omega_{k},g_{k},\Gamma_{k})$ and $w_{n}(\Omega_{k},g_{k},\Gamma_{k},c_{k})$ from the effective environment. This is in general a highly nontrivial inversion problem involving an underdetermined, non-convex system of nonlinear equations of mixed degrees, and can be hard to solve: there is a trade-off between this complexity and the accuracy of the initial fit, with an optimum at small numbers ($N\leqslant 5$ in all our applications) of interacting oscillators. Neither existence nor uniqueness of solutions are guaranteed for this inversion problem and physical requirements such as positivity of the rates $\Gamma_{n}$ need to be taken into account as well, so it is typically necessary to minimize some distance between $C^{R}(t)$ and $C^{E}(t)$ instead of exactly matching the best fit $\tilde{C}^{E}(t)$; this change in the correlation function is the only error introduced into the problem by the use of an effective environment.

In some cases with $N\leqslant 3$, it is possible to invert the equations exactly and obtain valid effective bath parameters; we have listed a few explicit solutions in Appendix B and will put some of them to use in our example applications. For general $N$, we devised a variational recipe to carry out our Transformation to Surrogate Oscillators (TSO) systematically. This is described in detail in Appendix A. The whole procedure can be summarized as follows:

• •

Fit $C^{E}(t)$ with $N$ complex exponentials $e^{\tilde{\lambda}_{n}t}$ with complex coefficients $\tilde{w}_{n}$ such that $\sum^{N}_{n=1}\tilde{w}_{n}>0$.

• •

Sample random points in a suitably sized open set $(0,g_{\mathrm{max}})^{N-1}$, to be used as coupling constants.

• •

Substitute each $(N-1)$-tuple $(g_{1},\dots,g_{N-1})$ into the equations relating the complex eigenvalues $\lambda_{n}$ to the $\tilde{\lambda}_{n}$ and solve: this will give rates $\Gamma_{n}$ and frequencies $\Omega_{n}$ such that the eigenvalues match; accept only solutions with all $\Gamma_{n}>0$.

• •

Compute the left and right eigenvectors of the matrices $M$ corresponding to each solution found, plug them into Eq. (21) and minimize a distance (e.g. the Manhattan distance $d_{\mathrm{Man}}(\mathbf{w},\mathbf{\tilde{w}})\coloneqq\sum^{N}_{n=1}|w_{n}-\tilde{w}_{n}|$) between $\mathbf{w}$ and $\mathbf{\tilde{w}}$ by varying the $c_{n}$.

• •

Assess overall accuracy of the solutions found and rank the corresponding $C^{R}(t)$ according to a meaningful figure of merit, such as the integral $\int^{t}_{0}\mathrm{d}t^{\prime}\int^{t^{\prime}}_{0}\mathrm{d}t^{\prime\prime}|C^{R}(t^{\prime}-t^{\prime\prime})-C^{E}(t^{\prime}-t^{\prime\prime})|$ from Ref. (SpinBosonBounds, ).

• •

If the accuracy of all effective correlation functions obtained is deemed insufficient, repeat with one more mode.

Effective environments obtained through this procedure can then be used to simulate the reduced dynamics of any model in which the interaction with the bath is mediated by the same bath operator $F_{R}$: the TSO is carried out once and for all irrespective of the system coupled to the environment given, and the effective environment can be used in any problem involving the same correlation function $C^{E}(t)$. We also wish to remark that for composite systems with multiple local environments, the procedure applies to each independent correlation function individually and yields local effective environments to be coupled to the corresponding parts of the system in the same way as the original ones, with no further complications arising: this feature is sketched in Fig. 3 and will be demonstrated in Section V.

By the same token, complex correlation functions requiring many exponentials for an accurate fit can be treated by breaking down the effective environment into smaller clusters of interacting modes, with each cluster accounting for a different component of $C^{E}(t)$—or equivalently, the underlying spectral density $J(\omega)$ of the unitary environment—as shown in Fig. 3. Note that decoupling all oscillators from each other, i.e. taking all $g_{n}=0$ (which corresponds to requiring all $\tilde{w}_{n}$ to be real and positive in Eq. (22)), one recovers the noninteracting pseudomodes of Ref. (Garraway_pseudomodes, ) as a limiting case.

III.3 Working example: Ohmic spectral density

To better illustrate the technique explained in the previous subsection, let us now demonstrate how our transformation works with an explicit example.

Consider an arbitrary quantum system coupled to an infinite environment in thermal equilibrium through the position operator of each oscillator (for a more succinct notation, we will leave the tensor products implicit and use natural units $\hbar=1$, $k_{B}=1$ from now on):

[TABLE]

This type of coupling for microscopic models is one of the most common in the OQS literature (BreuerPetruccione, ; Weiss, ; CaldeiraLeggettModel, ; Leggett_SpinBoson, ; DeVegaAlonso_NonMarkovianity, ). For a thermal initial state at inverse temperature $\beta=1/T$, the correlation function of the interaction operator $G_{E}=\int^{\infty}_{0}\mathrm{d}\omega\,g(\omega)(a_{\omega}+a^{\dagger}_{\omega})$ is

[TABLE]

where the spectral density $J(\omega)$ is related to the frequency-dependent coupling strength $g(\omega)$ through

[TABLE]

and typically given as a starting point for studying the problem. Spectral densities are real and positive by definition, and are often categorized according to the power of $\omega$ best approximating their behavior near the origin, where they are always zero; a $J(\omega)\propto\omega^{s}$ is called Ohmic if $s=1$, and super-(sub-)Ohmic if $s>1$ ($s<1$). The spectral density and temperature uniquely determine $C^{E}_{\beta}(t)$ and, consequently, the effect of the environment on the system.

Note that for unitary environments the correlation function is Hermitian in time, i.e. its real part is even and its imaginary part is odd, as can be seen clearly from Eq. (24). This implies that its Fourier transform

[TABLE]

where $\theta(\omega)$ is the Heaviside step function, is always real; at temperature $T=0$, it is just $2J(\omega)\theta(\omega)$. In fact, $C^{E}_{\beta}(\omega)/2$ may itself be regarded as a spectral density defined over a new environment, which comprises both positive- and negative-frequency modes and gives the correlation function $C^{E}_{\beta}(t)$ if initialized in the vacuum state (MayKuehn, ): this allows one to effectively rephrase arbitrary-temperature OQS problems as zero-temperature ones if it is convenient to do so, a possibility exploited by thermofield-based and other numerical methods (DeVegaBanuls_Thermofield, ; DiosiGisinStrunz_NonMarkovianity, ; Ritschel_AbsorptionSpectra, ; TamascelliSmirne_ThermalizedTEDOPA, ).

For non-unitary environments, in which time evolution is not an invertible map, correlation functions $C^{R}(t)$ are only defined at positive times; we extend the definition to negative times by imposing the same symmetry

[TABLE]

in order to be able to compare exact and effective correlation functions in the frequency domain instead of inspecting their real and imaginary parts separately.

Consider an Ohmic $J(\omega)$ with an exponential cutoff:

[TABLE]

Ohmic spectral densities define a very important class of environments entering the study of many systems, such as a particle undergoing quantum Brownian motion, or microscopic models leading to a Lindblad equation for a harmonic oscillator in a weakly coupled high-temperature environment (BreuerPetruccione, ; CaldeiraLeggettModel, ; FordLewisOConnell_DampedOscillator, ). The thermal correlation function $C^{E}_{\beta}(t)$ corresponding to the spectral density defined in Eq. (26) can be determined analytically as

[TABLE]

where $\psi^{\prime}(z)\coloneqq\frac{1}{\Gamma(z)}\frac{\mathrm{d}\Gamma(z)}{\mathrm{d}z}$ is the polygamma function of order one.

Performing our TSO on this correlation function at temperature $T=\Omega_{c}$ according to the recipe described in the previous subsection, for $N=4$ we determined the parameters given in Table 1; Fig. 3 shows the Fourier-transformed effective correlation function $C^{R}_{\beta}(\omega)$ obtained using these parameters and the target $C^{E}_{\beta}(\omega)$ for comparison. As can be seen from the plot, four interacting oscillators were enough to obtain a very accurate $C^{R}_{\beta}(t)$, with a peak in the error around $\omega=0$ reaching about 2% of the function value (see the inset of Fig. 3). This error affects the correlation function at very long times compared to its decay time, so we expect it to have a minor impact on the transient reduced dynamics of the system and become potentially more important at very long times. In all our tests, a small region around the origin was consistently found to be the part of the frequency domain where a general $C^{E}_{\beta}(\omega)$ is hardest to match: this is because any $C^{R}_{\beta}(\omega)$ is analytical around zero by construction, whereas $C^{E}_{\beta}(\omega)$ has discontinuous derivatives, as can be checked from Eq. (24). We stress again that the nonzero temperature is encoded in the effective parameters and not in the initial state, allowing us to treat very different temperature regimes at comparable costs, as will be made clearer in the next sections.

IV A test case: the spin-boson model

We now turn to the second part of our approach: computing the reduced dynamics of a system by coupling it to the effective environment and solving the relevant Lindblad master equation (8).

In order to demonstrate and quantitatively validate the method, we will show here the results we obtained for a system for which an analytical solution is known: the purely dephasing spin-boson model (BreuerPetruccione, ; Leggett_SpinBoson, ; CaldeiraLeggett_PureDephasing, ). The Hamiltonian for this system is

[TABLE]

and we consider again the Ohmic spectral density defined in Eq. (26). In this model, the system and interaction Hamiltonians commute and are both diagonal in the system basis, so the populations $p_{0}\coloneqq\rho_{00}$ and $p_{1}\coloneqq\rho_{11}$ are conserved by the evolution. Any coherence in this basis present in the initial state, on the other hand, is erased according to the law (see Ref. (BreuerPetruccione, ) for a derivation)

[TABLE]

with

[TABLE]

Using the cutoff frequency of the environment $\Omega_{c}$ as our energy scale, we set the parameter values $\omega_{0}=4\Omega_{c}$ and $k=1$, corresponding to a strong-coupling regime. Comparable coupling strengths appear e.g. in the study of superconducting quantum transmission lines (Peropadre_UltrastrongCoupling, ). The system is initialized in the pure state $\rho_{0S}=\ket{+}\!\bra{+}$, with $\ket{+}\coloneqq\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$ in terms of the eigenstates of $\sigma_{z}$, and we simulated the reduced dynamics at three different temperatures $T=0$, $T=\Omega_{c}$ and $T=\frac{5}{2}\Omega_{c}$. Recall that the effective bath is always at zero temperature; different temperatures of the original environment require different surrogate baths. We found accurate effective correlation functions with $N=4$ for the first two cases, and with $N=5$ for the high-temperature regime; the parameters are given in Appendix C, and the errors of the two correlation functions at nonzero temperatures are similar to the zero-temperature case already discussed.

IV.1 Results, accuracy and performance

We solved the effective Lindblad equations for all three cases using the QME integrator provided by the Python OQS package QuTiP (QuTiP1, ; QuTiP2, ), which implements a twelfth-order Adams-Moulton discrete integration algorithm.

From the results shown in Fig. 4, we see that our simulations with effective correlation functions give quantitatively good results for the coherence $\rho_{01}(t)$ at all times and temperatures (the populations $p_{0}$ and $p_{1}$ are both equal to $1/2$ throughout the evolution), as the overlap between the numerical (solid lines) and exact (dashed lines) solutions shows. The pure quantum decoherence at $T=0$ induces an algebraic decay asymptotically proportional to $t^{-k^{2}}$, while at $T>0$ the damping becomes exponential; a stronger effective coupling regime, which is determined both by $k$ and the strength of thermal effects, induces faster relaxation in the system dynamics.

The plots in the insets show the error figure

[TABLE]

which is identical for $f=\Re[\rho_{01}]$ (solid line) and $f=\Im[\rho_{01}]$ (dashed line). This is a better estimator for the accuracy of the simulations than e.g. the absolute difference $|f(t)-f_{\mathrm{Num}}(t)|$ because it removes the bias coming from changes in the relaxation time due to temperature, allowing us to compare all regimes on an equal footing. The error, as measured by $E_{f}(t)$, remains of the order of a few percent until the system has almost reached equilibrium and is comparable for the three regimes probed, mirroring the similar relative errors we had in all three effective correlation functions. The latter observation can be understood as follows: as higher temperatures or larger coupling constants increase the effects of the bath on the system, the error being carried from the correlation function into the reduced dynamics is magnified accordingly; on the other hand, these stronger effective regimes shorten the relaxation time of the system, so the cumulative effect of the error over time is not as severe as when the coupling is weaker or the temperature is lower.

From these results, we conclude that the method is quite reliable and stable provided that the effective correlation functions used are reasonably accurate, and that this accuracy does not command significantly greater effort or complexity in the TSO at higher temperatures and is completely independent of the system and the coupling strength. Furthermore, it is worth noting that any method based on approximating the environment alters its correlation function and is therefore prone to the same kind of error as ours, but we use a rigorously motivated and physically meaningful quantifier to optimize our correlation functions and keep it under control.

The computational cost of the simulations depends on the local dimension at which each effective mode is truncated and on the spread between the total evolution time and any faster timescales in the problem at hand, though the memory requirements scale faster with the complexity of the problem than the computation times do; to obtain our converged results for this system, which required $\sim 4$ levels for most oscillators and a maximum of $7$ for one mode in each simulation, all running times were below 10 minutes on a laptop. This cost grows rapidly with the number of effective modes, the local dimensions needed for convergence and the size of the system itself; on the other hand, temperature and coupling strength have a limited impact on these factors: a very strong coupling or high temperature will require higher local dimensions but also cause very rapid relaxation to equilibrium, making long simulation times unnecessary. Moreover, when the coupling is stronger and more levels are needed for convergence, this typically affects one particular mode much more than the others, leading to an effective polynomial rather than exponential scaling in the coupling strength and temperature.

V Physically relevant applications

Few- and many-body systems non-perturbatively coupled to non-Markovian environments with structured spectral densities are ubiquitous in many fields ranging from biological physics and chemistry (HuelgaPlenio_QBio, ; Pelzer_Transport, ; DeSio_OPV, ) to condensed matter (RibeiroVieira_Transport, ; Haase_Metrology, ), thermodynamics (UzdinLevyKosloff_QHE, ; MitchisonPlenio_NonEquilibrium, ), nanomaterial science and sensing (JelezkoPlenio_NV, ) and quantum metrology (ChinHuelgaPlenio_Metrology, ; SmirneKolodynski_Metrology, ; HaaseSmirneKolodynski_Metrology, ), and have prompted much research in theoretical modeling and numerical simulation methods for general OQS. In this section, we will demonstrate how our approach may be used to solve and make predictions on models at the forefront of current research, by presenting the results we obtained in two different applications.

In the first of the two examples, we will show some results for experimentally measurable optical properties in a model inspired by studies on coherent charge and energy transfer in biological molecular aggregates (Scholes_QBioNature, ). After that, we will consider a polymer model of the type relevant in research on organic photovoltaic materials (Clark_OPV, ; Tamura_OPV, ), and report simulation results for excitation transport dynamics in such a system under the assumptions of strongly coupled, non-Markovian local environments interacting with each monomer. The former example demonstrates the use of our simulation technique to gain physical insight by direct comparison of the results with observable data, while the latter gives an idea of its potential in terms of performance by addressing a problem beyond the reach of other current state-of-the-art methods.

V.1 Optical spectra in molecular aggregates

We considered a simple dimer model with system parameters in the range of those found in biomolecular aggregates participating in excitation energy transfer (Plenio_QBio, ; PlenioAlmeidaHuelga_Dimer, ), coupled to an environment with a realistic spectral density derived from models common in the literature (AdolphsRenger, ). We first simulated the reduced dynamics of this system at liquid nitrogen temperature ($T=77\,\mathrm{K}$), comparing the results with a simulation done using the numerically exact TEDOPA (PriorPlenio_TEDOPA, ; TamascelliSmirne_ThermalizedTEDOPA, ; Tamascelli_RSVD, ; Kohn_RSVD, ), and then computed its absorption spectrum for the same temperature as well as for $T=0\,\mathrm{K}$ and $T=300\,\mathrm{K}$. Two different simulation techniques were used to integrate the effective Lindblad equation for the dynamics and the absorption spectra, and the spectra were calculated for two different environmental spectral densities and compared in order to identify the optical signatures setting them apart: in particular, we sought to determine the differences between spectra obtained in the presence or absence of a strongly coupled vibrational mode in addition to a broad background spectral density.

V.1.1 Details of the model and reduced dynamics

Following Ref. (PlenioAlmeidaHuelga_Dimer, ), we considered a free dimer Hamiltonian

[TABLE]

where the two monomers have on-site energies $E_{1}=E_{g}+12\,328\,\mathrm{cm}^{-1}$ and $E_{2}=E_{g}+12\,472\,\mathrm{cm}^{-1}$ and interact through a hopping coefficient $J=70.7\,\mathrm{cm}^{-1}$, and $\ket{E_{12}}$ is the state with both monomers excited. We only considered the ground state $\ket{g}$ and the single-excitation manifold spanned by the states $\ket{E_{n}}$, ignoring the doubly excited state $\ket{E_{12}}$ since its contribution is typically negligible both in excitation transport phenomenology and most absorption experiments (MayKuehn, ). Then, setting $E_{g}=0$ as our reference energy, we are left with an effective two-level Hamiltonian for the single-excitation manifold

[TABLE]

whose eigenstates $\ket{\varepsilon_{1,2}}$ have an energy gap of $\Delta=201.8\,\mathrm{cm}^{-1}$, with the ground state dynamically decoupled and only contributing to expectation values or correlation functions of operators explicitly dependent on it.

The local excited states $\ket{E_{n}}$ interact with separate environments, which account for the molecular vibrations (both within the system and in the protein scaffold around it) and the presence of a solvent. We model these degrees of freedom by coupling the monomers to independent thermal baths with the same spectral density and temperature; the physical model is sketched in Fig. 5 (a).

We first studied the problem for a spectral density consisting of two contributions: a broad background noise spectrum in the super-Ohmic form first introduced by Adolphs and Renger (AdolphsRenger, )

[TABLE]

where the two cutoff frequencies are $(\Omega_{\mathrm{AR}1},\Omega_{\mathrm{AR}2})=(0.557,1.936)\,\mathrm{cm}^{-1}$ and the weights of the two terms are $(\rho_{1},\rho_{2})=\frac{288}{5}(\frac{8}{13},\frac{5}{13})$, and a strongly coupled vibrational mode represented by adding an antisymmetrized Lorentzian peak

[TABLE]

For this sharp spectral feature, we set $\Omega=227.5\,\mathrm{cm}^{-1}$, slightly above resonance with the system, a width $\Gamma=20\,\mathrm{cm}^{-1}$ corresponding to a decay time $(\Gamma/2)^{-1}\sim 0.5\,\mathrm{ps}$, and a Huang–Rhys factor $S=0.0379$ placing it in a moderate-coupling regime with the system. The reorganization energies corresponding to the background and the full environment are

[TABLE]

In order to compute the reduced dynamics of the system in the single-excitation subspace, the total Hamiltonian of our problem

[TABLE]

can be rewritten in terms of the ‘common-mode’ and ‘relative’ creation and annihilation operators parametrized by a single frequency $A^{(\dagger)}_{\omega}=\frac{a^{(\dagger)}_{\omega_{1}}+a^{(\dagger)}_{\omega_{2}}}{\sqrt{2}}$ and $a^{(\dagger)}_{\omega}=\frac{a^{(\dagger)}_{\omega_{1}}-a^{(\dagger)}_{\omega_{2}}}{\sqrt{2}}$:

[TABLE]

The common-mode environment only interacts with the single-excitation subspace through the last term, which is proportional to the identity in that subspace. Therefore, it can be ignored in any calculation not involving the ground state: for such applications, the Hamiltonian then reduces to

[TABLE]

in terms of the relative modes only, and the dynamics factorizes between the two subspaces unless coherences between them are present in the initial state. A sketch of the model after this rearrangement of the environmental modes and the TSO is given in Fig. 5 (b).

We computed the reduced dynamics in the single-excitation subspace for an initial coherent superposition of energy eigenstates $\rho_{0S}=\ket{\tilde{+}}\!\bra{\tilde{+}}$, where

[TABLE]

for the spectral density $J(\omega)=J_{\mathrm{AR}}(\omega)+J_{\mathrm{AL}}(\Omega,\Gamma,S;\omega)$ considered. To this end, we determined effective parameters corresponding to the two terms of $J(\omega)$ and temperatures $T=0\,\mathrm{K}$, $T=77\,\mathrm{K}$ ($53.5\,\mathrm{cm}^{-1}$) and $T=300\,\mathrm{K}$ ($208.5\,\mathrm{cm}^{-1}$), performing the TSO separately on the Adolphs–Renger background, Eq. (33), and the antisymmetrized Lorentzian peak, Eq. (34). This corresponds to assigning a separate effective environment to each additive part of the spectral density $J(\omega)$ and can be a convenient strategy to break down structured spectra, as mentioned in an earlier section and shown in Fig. 3. The Adolphs–Renger correlation function required $N=4$ oscillators at all three temperatures, and the Lorentzian mode was replaced by one effective oscillator at $T=0$ and two interacting ones at $T>0$ using the exact methods for $N=1,2$ described in Appendix B. All parameters of the effective environments are given in Appendix C. The environmental correlation function at $T=77\,\mathrm{K}$, the temperature for which we computed the dynamics, is plotted along with the effective correlation function from the TSO in Fig. 7. The other temperatures will be considered in the calculation of absorption spectra for the model dimer.

Since the amount of memory required for a direct integration of the effective Lindblad equation would become too large for the system coupled to six effective modes with the local dimensions needed for convergence, we carried out the simulations using the quantum jump or Monte Carlo Wave Function (MCWF) method for pure states (DumZollerRitsch_MCWF, ; DalibardCastinMolmer_MCWF, ; PlenioKnight_MCWF, ) instead (the memory cost of MCWF scales linearly with the total Hilbert space dimension $\mathcal{N}$ for sparse Lindblad superoperators such as ours, while a master equation integrator requires at least $O(\mathcal{N}^{2})$). The simulation was performed using another QuTiP code, since the package also provides MCWF routines. Our averages converged after as few as 1000 trajectories (this is due to the quantum jumps in the evolution directly affecting only the modes but not the system, since the latter has no Lindblad damping of its own, and thus partly canceling in the trace); we computed twice as many trajectories as a check but found no visible differences. The results of our simulation are shown in Fig. 7 along with those obtained by using TEDOPA: again, the accuracy of our effective correlation function—with errors of the order of 1% as in the previous section—translates to a satisfactory result for the reduced dynamics throughout the time window considered, which is almost enough for the system to reach equilibrium (no comparison was possible for times longer that about $1.3\,\mathrm{ps}$ due to the rapidly increasing cost of the TEDOPA simulation at later times). The numerical cost is also remarkably low: for converged local dimensions, the simulation required under $200\,\mathrm{MB}$ of memory per thread and could therefore have been carried out on a desktop or laptop computer. To achieve higher parallelization of the work, however, we used the JUSTUS cluster at Ulm University: on a 16-core cluster node, the reduced dynamics up to $t=2\,\mathrm{ps}$ took 22 minutes to compute and the scaling is linear in the total simulation time. For comparison, TEDOPA took around 60 minutes to reach $t=1.3\,\mathrm{ps}$ on the same hardware and started scaling superlinearly in the simulation time at around that point.

V.1.2 Absorption spectra

Absorption experiments probe the linear response of the system; light from a laser source can be described as interacting with the local dipole moment operators $\vec{\mu}_{n}\coloneqq\vec{d}_{n}\ket{E_{n}}\!\bra{g}$, where $\vec{d}_{n}$ is the classical dipole moment of the $n$-th site, in a perturbative manner (Carmichael, ; MayKuehn, ). Then the spectrum is obtained from the one-sided Fourier transform of the correlation function of the total dipole operator $\vec{\mu}\coloneqq\sum^{2}_{n=1}\vec{\mu}_{n}$ over the initial stationary state

[TABLE]

where the bath is in a thermal state at inverse temperature $\beta$ and the system is in the electronic ground state, which does not couple to the environment, since excited-state populations at equilibrium are negligible due to the very low intensity of the laser in such a setup (Mukamel, ; MayKuehn, ).

Specifically, the correlation function of interest is given by the scalar product of the dipole operator $\vec{\mu}$, applied at times $t_{0}=0$ and $t$: in terms of the overall unitary evolution, one has

[TABLE]

Note that this is formally a two-time object: we can compute it using an effective environment because the first operator acts on the system at equilibrium, so the hypotheses of Theorem 1 are not violated. The unitary dynamics acts on $\left(\vec{d}_{1}\ket{E_{1}}+\vec{d}_{2}\ket{E_{2}}\right)\!\bra{g}\rho_{0\mathrm{Abs}}$, which is still a factorized object with the environment in a thermal state, so the equivalence with a suitable effective Lindblad dynamics remains well defined; however, note that this time the common-mode part of the total environment does not decouple from the problem and one needs to simulate the system along with both local baths, as pictured in Fig. 5 (c).

We set the ansatz $\vec{d}_{1}=\vec{d}_{2}=\vec{d}$ for the geometry of the dimer in order to simplify the form of the dipole correlation function. Expressed in units of $|\vec{d}|^{2}$, $C_{\mu}(t)$ becomes

[TABLE]

The absorption spectrum is then given by

[TABLE]

In order to compare the effect on the absorption spectrum of a strongly coupled, underdamped vibrational mode in the environment, we will now consider two spectral densities: $J_{0}(\omega)\coloneqq J_{\mathrm{AR}}(\omega)$ and $J_{1}(\omega)\coloneqq J_{\mathrm{AR}}(\omega)+J_{\mathrm{AL}}(\Omega,\Gamma,S;\omega)$, with $\Omega=200\,\mathrm{cm}^{-1}$, $\Gamma=10\,\mathrm{cm}^{-1}$ and $S=0.25$. A plot of the spectral densities is given in Fig. 8: as the figure shows, the contribution of the underdamped peak is much stronger in this new setup.

Integrating the effective Lindblad equation with the initial pseudo-state $\tilde{\rho}_{0}\coloneqq\left(\ket{E_{1}}+\ket{E_{2}}\right)\!\bra{g}\rho_{0\mathrm{Abs}}$, one obtains the dipole correlation from Eq. (40) as

[TABLE]

where $\mathcal{L}$ is the Lindblad superoperator given by the TSO with both local environments included.

Computing the dipole correlation function is much more challenging than simulating the reduced dynamics in the single-excitation subspace, because this time both sets of surrogate modes need to be explicitly accounted for and the local dimensions are quite high, as shown in the relevant parameter tables in Appendix C. In order to keep the total Hilbert space dimension manageable, we employed a variation of a recently published tensor network–based technique called Dissipation-Assisted Matrix Product Factorization (DAMPF) (Somoza_DAMPF, ) to simulate the dimer. DAMPF, which was originally developed using non-interacting pseudomodes, is extremely efficient for vibronic aggregates in the single-excitation manifold, making it an ideal candidate for scaled-up simulations of systems involving many sites with local surrogate-oscillator baths, as will be shown in the next subsection.

For each temperature, we simulated the dimer until $C_{\mu}(t)$ had decayed to values small enough for the limit in Eq. (41) to be approximately satisfied (the initial value in our units is $C_{\mu}(0)=2$, as can be seen from Eq. (42)): the final times and corresponding absolute values of $C_{\mu}(t)$ reached in our simulations are reported in Table 2, and the resulting absorption spectra—obtained via a discrete Fourier transform and centered around the midpoint frequency $12\,400\,\mathrm{cm}^{-1}$ of the single-excitation subspace—are shown in Fig. 9.

The spectra show the expected features: the result for $J_{0}(\omega)$ displays absorption lines corresponding to the single-excitation eigenstates $\ket{\varepsilon_{1,2}}$ of $H_{\mathrm{1ex}}$ and appearing at the corresponding energy values redshifted by the bath reorganization energy; the line corresponding to the higher eigenstate is broadened due to the decay channels of that state, which couples to the environment and the lower excited state, whereas the latter gives a very narrow zero-temperature peak since it is not coupled to any lower-lying state it could decay to. At higher temperatures, the contribution from the upper level becomes larger than the lower one, but the energies associated with the—now markedly broadened—spectral lines no longer represent energy eigenstates of the system, since the dressed system-environment energy eigenbasis is very different from a tensor product basis in this regime, as hinted at by the fact that the environmental reorganization energy corresponding to the thermalized spectral density $C^{E}_{\beta}(\omega)/2$ is comparable to $J$. At room temperature, hardly any structure is discernible but for the fact that the spectrum rises slowly and somewhat irregularly to the left of the maximum.

Adding the strong peak to the spectral density, the spectra are shifted to the left by the added reorganization energy, and the expected new spectral lines associated with excitations of both the dimer and the coupled vibrational mode appear. At lower temperatures, higher sidebands are also visible as small bumps to the right of the main spectral curves; they are washed out by the strong broadening at room temperature.

It should be noted that at temperatures up to $T=77\,\mathrm{K}$ the timescale at which the reduced dynamics of the system reaches a steady state is of the order of about one picosecond for both $J_{0}(\omega)$ and $J_{1}(\omega)$ (most of the dissipation is due to the broad Adolphs–Renger background, since the Lorentzian mode has a long lifetime); the decay times of the dipole correlation functions, on the other hand, were found to be significantly longer. In order for $C_{\mu}(t)$ to reach values small enough to avoid visible spurious effects from an incomplete decay in the Fourier transform, some of the simulations had to run up to times of order $t=10\,\mathrm{ps}$ (see again Table 2). Such long-time simulations are only possible with methods whose cost scales slowly, e.g. linearly, in the simulation time, such as DAMPF (or, for smaller Hilbert spaces, MCWF or even direct integration routines). Convergence in DAMPF is achieved for sufficiently high local dimensions as well as bond dimensions $\chi$ (we refer the reader to Ref. (Schollwock_tDMRG, ) for more details on parameters in the tensor-network setup); we found that both needed to be quite high at nonzero temperatures, and used a time-adaptive bond dimension for the $T=300\,\mathrm{K}$ case to optimize the effective Hilbert space dimension throughout those simulations in order to save time. The local dimensions of each mode are given in the relevant parameter tables in Appendix C, and the bond dimensions are shown in Table 2 for those simulations in which they were kept fixed.

V.2 Excitation transport in organic polymers

In this subsection, we will give a more concrete demonstration of the full potential of surrogate environments for physically sound and numerically efficient simulation of systems. To this end, we will apply our method to a problem involving an organic polymer modeled as a chain consisting of many sites with realistic local environments strongly coupled to each of them.

Organic polymers have been gaining growing attention from the condensed- matter, OQS and many-body-physics communities due to their considerable technological potential, e.g. in devising novel photovoltaic and other electronic components (Clarke_OPV, ; Proctor_OPV, ). Such systems are often modeled in the same tight-binding approximation used for photosynthetic complexes in biological physics, and simulating charge transfer or separation processes in chains of organic monomers interacting with local non-Markovian environments is a notoriously challenging task even with state-of-the-art techniques such as HEOM, as mentioned previously (Yamagata_QuantumWires, ; Tempelaar_Coherence, ; Spano_Aggregates, ; Hestand_Aggregates, ).

Typical treatments of such organic systems often employ strong coarse-graining of the environmental spectral features (Blau_dimer, ; Chenel_OPV, ), in order to save computational resources for the simulation of a system which may consist of a large number of sites. We will now show how our surrogate environments can be used to calculate the reduced dynamics of extended vibronic systems consisting of multiple sites, with each site coupled to a realistic thermal bath comprising both a sharp mode and an Ohmic background. Interactions, both among sites and between each site and its local environment, are strong, with the baths characterized by high reorganization energies, and we will consider the system at both zero and room temperature.

We considered a homogeneous polymer Hamiltonian of the form

[TABLE]

where $K$ is the number of sites, $J$ is the site-site coupling and the on-site energies are assumed equal and set to zero. The ground state is disregarded since it does not couple to the single-excitation subspace we are working in, and we set $K=10$ and $J=200\,\mathrm{cm}^{-1}$.

Each site couples to a local thermal bath in the same way as in Eq. (35), and the spectral density of the baths is a sum of an Ohmic background of the form Eq. (26), with cutoff frequency $\Omega_{c}=200\,\mathrm{cm}^{-1}$ and rescaled by an overall factor $\kappa=0.25$, and an underdamped peak at $\Omega=1\,000\,\mathrm{cm}^{-1}$ with $\Gamma=20\,\mathrm{cm}^{-1}$ and Huang–Rhys factor $S=0.25$. The total reorganization energy is

[TABLE]

a very high value.

We simulated the evolution of this system up to $t=1.25\,\mathrm{ps}$ from an initial state $\rho_{0S}=\ket{E_{1}}\!\bra{E_{1}}$ at $T=0\,\mathrm{K}$ and at $T=\Omega_{c}=288\,\mathrm{K}$, again separating the background and the peak in the TSO. For the Ohmic spectral density, we used the surrogate environments already introduced in Sections III and IV. The results are shown in Figs. 11 and 11, respectively. The local dimensions for the Ohmic background needed to be higher for this simulation than for the spin-boson case discussed in Section IV (the population dynamics in the surrogate oscillators depends both on the coupling strength and on the internal dynamics of the system they are interacting with in any given problem), and we saw no relevant changes in the results for $\chi>9$ at zero temperature and $\chi>12$ at room temperature. The zero-temperature simulation took a few hours and the room-temperature one was completed over the course of several days on a 16-core node of the JUSTUS cluster at Ulm University.

Simulations such as these on a conventional architecture (i.e. one not boosted by the use of graphical processing units, which would further enhance the numerical efficiency of both other schemes and our own) are beyond the reach of any simulation technique we know of. The strength of the coupling, especially to the high-frequency mode, would be critical for a HEOM treatment with as few as two sites (Somoza_DAMPF, ), and the system size rules out TEDOPA, QUAPI or any other method for non-Markovian open systems regardless of that property of the environment. Regarding DAMPF, which overcomes the problem of ever-growing bond dimensions in tensor network–based methods, it can also be used with independent auxiliary oscillators. However, this comes at the price of using a much greater number of modes for the same accuracy, again driving up the simulation cost. With the interacting modes given by the TSO, we can be certain that our results are closer to the reduced dynamics of the unitary model than anything that can be done with the same number of independent pseudomodes, due to the far smaller correlation function error our coupled modes entail. Therefore, this combination of accuracy and numerical efficiency would not be possible otherwise.

VI Discussion

After introducing and demonstrating our new simulation method, let us now recapitulate its main theoretical and technical points, discussing its strengths, limitations and error sources in order to give a clear and concise summary of its current state and possible future improvements.

VI.1 Theoretical basis and general remarks

Our method is part of a category of hybrid approaches based on rephrasing microscopic OQS models as effective Markovian problems, in which the memory of the environment is accounted for by an ad hoc auxiliary system. Although this divide-and-conquer strategy between Markovian and non-Markovian effects is a shared feature of several existing methods, the flexibility and quantitative control allowed for by the rigorous theoretical groundwork underlying our construction (TSO_Theorem, ; SpinBosonBounds, ) are, as far as we know, unprecedented for an approach of this type.

The transformation procedure we described in Section III exploits the generality of a broad, physically well-defined class of effective environments to tailor them in a systematic way to fit the microscopic ones given: isolating the correlation function as the single property of the environment which needs to be replicated as accurately as possible, we take advantage of the added versatility from using interacting effective modes to make this fitting procedure more accurate while keeping the number of effective degrees of freedom lower than would be possible in chain- or star-configuration schemes. It should be noted that we took one particular ansatz for the effective environment because we found it the most convenient for our needs, but many other choices (non-hopping linear couplings, interactions beyond nearest neighbors, different damping and initial stationary state, etc.) are possible.

Another relevant feature is that the system on which the environment acts does not enter in this part of the procedure at all. Therefore, once the effective parameters corresponding to a given unitary bath are determined, they can be used in all problems featuring that particular bath, as we showed in Section V. This makes determining the surrogate environment a one-off task, which can be very convenient in any field in which standard spectral densities appear in many different situations.

The second part of the method is the simulation of the system coupled to the effective environment. Here any of the analytical or numerical techniques for Lindblad master equations already developed in the literature can be used; the system part of the problem is completely unrestricted, so different strategies can be adopted depending on the problem at hand and available computational resources. We demonstrated the method on two-level systems interacting with small sets of up to six effective modes, for which simple and clear solution methods like direct integration of the master equation or MCWF are still suitable, and on more complex systems such as chains of monomers in the tight-binding approximation, with each local site interacting with its own surrogate bath, for which an integration scheme compressing the total Hilbert space dimension is necessary.

In summary, we have shown that although the mathematical question of how to generalize the use of independent damped oscillators as effective environments is highly nontrivial, finding ways to do so can be extremely beneficial to the modeling of non-Markovian OQS. Our recipe for the construction of few-body Gaussian environments with interacting surrogate modes proved a valuable technique to encode complex environmental effects in surprisingly compact effective models, with remarkable computational advantages.

VI.2 Numerical complexity and costs

As we have shown in the examples given in the preceding sections, our method is very versatile and applies in principle to any non-Markovian, non-perturbative OQS problem involving a Gaussian bosonic bath, at any temperature and coupled arbitrarily strongly to a system. Clearly, some problem classes and physical regimes are more suited than others to this type of treatment: we will now summarize the key elements determining the computational effort required for any given application.

Concerning the simulation of systems, the most important variable to look at is the Hilbert space dimension, which includes both the system size and the number and local dimensions of the surrogate modes; temperature and coupling strength affect the overall complexity indirectly, mainly by determining the minimum local dimensions needed to accurately simulate the modes. The number of effective oscillators and their parameters come from the TSO and depend on the spectral density and temperature of the unitary environment. The spectral density has a prominent role in determining the number of oscillators required; higher temperatures can contribute too but typically result in an effectively stronger coupling to the system instead: this does cause the oscillators to become more populated, making higher local dimensions necessary for the reduced dynamics to converge, but the added computational cost is usually less than that entailed by adding a new mode.

When simulating highly structured systems or environments, the Hilbert space dimensions involved are such that memory, rather than time, typically becomes the main computational concern. Hilbert space dimensions of order $10^{5}$ require memory of order $1\,\mathrm{GB}$ per thread with the MCWF implementation we used in the aforementioned simulation, which can be managed on desktop-level hardware or individual nodes of a cluster. Larger Hilbert spaces, such as those of the dimer and polymer problems we considered, must be compressed by suitable optimization techniques such as matrix product operator–based methods (Mascarenhas_MPO, ; CuiCiracBanuls_MPO, ; Somoza_DAMPF, ) or reduction to Krylov subspaces, a topic of current relevance in the study of large systems of numerical differential equations (Minchev_MatExp, ; VoSidje_Krylov, ; Tokman_KIOPS, ). We used the newly developed DAMPF technique because it could be easily adapted from its original form in order to accommodate coupled surrogate modes, and showed that the scaling of the simulations with the size and complexity of the system studied is quite favorable, in some cases outperforming any known method and thus attaining results hitherto out of reach.

In general, the cost of simulating a Lindblad dynamics scales linearly in the total evolution time for methods such as direct integration or MCWF, which allow for the full Hilbert space dimension to be fixed upfront. This makes them well suited for the study of long-time dynamics and relaxation to equilibrium. When the total effective Hilbert space of a problem is too large for any such technique, one needs to resort to time-adaptive truncation schemes, which can scale quite unfavorably in time. However, novel methods such as DAMPF exploit the damping in the simulated dynamics to bound the maximum effective dimensions they use, thus reducing these nonlinear additional costs to the point of recovering an approximately linear scaling which allowed us to simulate even a polymer up to arbitrarily long times. As ambiguous as performance assessments can become, depending on physical regimes and scales in the models studied, it seems quite clear nonetheless that there are situations in which using surrogate modes to reduce the number of effective degrees of freedom needed for accurate results is of paramount computational advantage.

As to determining the parameters of the surrogate modes in the first place, the inversion problem from the target correlation function is, in general, a mathematically difficult task. Our TSO algorithm uses a randomized parametrization as a variational method to reduce the number of variables in the problem and unlock a part of the solution, which is then fed back into the inversion problem to determine the values still missing; the solution found is the best possible for the random initial values given, and a minimization on the sample according to a suitable figure of merit is carried out a posteriori.

This rather involved procedure gives satisfactory results but scales poorly with the number of modes; for more than five interacting oscillators, it is already very expensive. This, however, is not a major setback for several reasons. First of all, complex environmental correlation functions typically originate from spectral densities comprising several simple terms, which can be addressed—and recycled for other problems if needed—individually, as demonstrated in our example applications; secondly, keeping the number of effective oscillators as low as possible is also a priority for simulation purposes and does not put significant constraints on accessible coupling or temperature regimes; finally, the cost of the TSO is not fixed but depends on the form chosen for the effective environment, so the complexity of our particular algorithm is not universal.

Possible future improvements to the variational algorithm could come from employing different methods such as simulated annealing or importance sampling in the parameter search or machine-learning techniques to minimize the distance between original and effective correlation functions with respect to the parameters; finding a way to work with the map from effective environment to correlation function in the direct rather than the reverse direction, if possible, would be a major simplification.

VI.3 Accuracy and error sources

To complete the discussion of our method, we must now turn to the sources of error affecting the reduced dynamics, and the control we have over them.

The most important error to be addressed is of physical origin and comes from the TSO. This is the error in the correlation function, and its impact on the reduced dynamics and operator expectation values at any time is rigorously bounded (SpinBosonBounds, ) (the paper focuses on the spin-boson model in particular, but similar bounds can be derived for other finite systems following the same steps).

Though under control, this error is worth a more careful analysis because it is actually a sum of two errors, one from a fundamental feature of our method, the other from a technical constraint.

The former source is the very form of any correlation function defined as in Eq. (10): it has been shown that no infinite, unitary thermal environment can have a correlation function of the form obtained via the quantum regression formula for a finite, Lindblad-damped auxiliary bath, because the fluctuation-dissipation theorem (CallenWeltonFDT, ; BreuerPetruccione, ), which holds for continuous, unitary thermal environments, is incompatible, strictly speaking, with the regression hypothesis (Talkner_NoQRT, ; FordOConnell_NoQRT, ). This is reflected quantitatively in the fact that no zero-temperature correlation function obtained through the latter is exactly zero on the whole negative frequency domain; however, the violation of the fluctuation-dissipation theorem can always be reduced by adding effective modes, until the unavoidable residual error is comparable to other errors in the model at hand (except in pathological cases, such as the weakly coupled spin-boson model with pure dephasing: this system is only sensitive to the limit of the spectral density at zero frequency, where the analytical differences between unitary and effective environments emerge most clearly, as discussed in Section III).

The second source contributing to the correlation function error are the constraints on the parameters in the master equation: for a fixed number of modes, not every linear combination of complex exponentials can be derived from a valid set of effective bath parameters via Eqs. (20) and (21) (for example, expressions obtained by setting at least one of the master equation rates $\Gamma_{n}$ to a negative value are out of physical scope). Therefore, the closest physically possible correlation function to the one given is generally not the best unconstrained fit with complex exponentials.

It should also be mentioned at this point that spectral densities of microscopic models are ultimately derived from experimental results in many applications of current interest (AdolphsRenger, ; Bennett_Dimer, ; Pelzer_Transport, ; JelezkoPlenio_NV, ), so any error in our TSO resulting in a correlation function still compatible with the data is immaterial in practice.

Finally, the last error source in our method is strictly numerical and comes from the integration schemes used to solve the Lindblad equation, which necessarily involve some truncation of the Hilbert space. This error is not under rigorous control, but requires method-dependent convergence checks like any other numerical solution technique.

VI.4 Impact

Finally, let us sum up the salient features of our simulation method and highlight its distinguishing qualities among existing schemes for general OQS.

First of all, we wish to emphasize that our aim in proposing this approach is to offer the level of accuracy and reliability of a fully microscopic simulation while retaining the benefits of working with two-tiered effective environments, particularly their simple mathematical structure and efficient numerical simulation.

Our scheme fills the gap between exact and simplified effective methods by providing auxiliary environments with a quantitatively certified link to the microscopic ones they stand in for, and enables very efficient simulation of nontrivial environments by keeping the number of modes much lower—thanks to the interactions among them—than any similar techniques we are aware of. For example, noninteracting pseudomodes are a special case of our surrogate oscillators, but we found that in order to reproduce an environment such as the one we discussed in our dimer example application in Section V, we would have needed at least 20 pseudomodes to attain the same accuracy given by our TSO with 6 oscillators: while each independent pseudomode contributes a term

[TABLE]

with a real and positive $w_{n}$—a Lorentzian, in the frequency domain—to the correlation function, our interacting effective modes contribute terms of the form

[TABLE]

thanks to the fact that the $w_{n}$ are complex, and these functions turn out to be far more flexible for fitting purposes. This difference is even more dramatic with Ohmic environments such as the one we considered in the polymer simulations, because the frequency-domain correlation function for an Ohmic bath at high temperature has a maximum at the origin: fitting such a shape with Lorentzians would result in a large number of underdamped low-frequency pseudomodes, which would become highly populated during the dynamics and require impractically high local dimensions. An accurate simulation of both our example systems would thus have been much more expensive using noninteracting pseudomodes.

We have also compared our simulations with calculations performed using microscopic methods and found that our accuracy is on par with numerically exact results, e.g. from TEDOPA, at least for the relatively simple application we used as a benchmark: long-time dynamics are much easier to compute by solving our effective Lindblad equation in all coupling and temperature regimes due to the nonlinear scaling of TEDOPA in the evolution time; on the other hand, spectral densities with complicated shapes requiring a large number of effective oscillators are hard with our method (though the solution of our effective Lindblad equation can be optimized by techniques such as DAMPF) while TEDOPA is much less sensitive to the shape of the spectral density. In our example application, a simple integration method running on a laptop performed better than TEDOPA at all temperatures for medium to long times (at zero temperature even for very short times) and comparably well for short times at nonzero temperature, for a moderately structured spectral density; the scaling in system size and complexity is very similar for the two schemes.

The HEOM method (Tanimura_HEOM, ; TanimuraKubo_HEOM, ) is more akin to our approach in spirit, since it is also based on exponential fitting of $C^{E}(t)$. Much like our number of surrogate modes, the number of exponentials needed for an adequate fit is one of the main factors determining complexity of HEOM simulations, the other being the tier at which the hierarchy needs to be truncated. This number should be the same for the two methods if one requires the same accuracy and uses efficient expansion techniques for the correlation function (Duan_EfficientHEOM, ; Duan_EfficientT0HEOM, ) (these overcome the well-known problem of the more traditional Matsubara-frequency expansion (MeierTannor_DecompositionHEOM, ; HuLuo_DecompositionHEOM, ), which at low temperatures needs a very large number of exponential terms in order to converge). HEOM scales with the complexity of the spectral density in a similar manner as our method and can likewise account implicitly for temperature through the approximate correlation function. Long evolution times are also not problematic for most regimes; however, they can be in the presence of certain environmental features, e.g. narrow peaks with high Huang–Rhys factors corresponding to strongly coupled environmental modes, which make the hierarchy of equations converge very slowly, significantly increasing the simulation cost. We have shown in the last of our example simulations that including the effect of such terms in our effective baths does not affect our simulation costs as dramatically as it does HEOM’s; in fact, the spectral density we considered poses a serious challenge for HEOM even with just two sites. Combined with the scaling of HEOM in the size of systems such as our dimer and polymer, this singles out at least one class of problems where both high- and low-temperature simulations make surrogate modes the most efficient if not the only viable option.

VII Conclusions and outlook

In this paper, we introduced a new non-perturbative approach for the description and simulation of arbitrary open quantum systems in Gaussian bosonic environments. The method is based on the use of networks of dissipative auxiliary oscillators as a means to account for nontrivial environmental effects, and puts no restrictions on temperature, non-Markovianity, system-environment coupling strength or system structure. It generalizes previously existing schemes employing independent fictitious modes, and we demonstrated that such a generalization is both sensible from a methodological point of view and extremely useful in terms of practical results.

We devised a systematic recipe to build effective environments of a very versatile class, tuning their parameters in order to capture the effects of any given unitary bath, using as few degrees of freedom as possible and with a clear measure of the error involved. This procedure is grounded in rigorous theoretical results, specifically the equivalence conditions between unitary and dissipative Gaussian baths proved in Ref. (TSO_Theorem, ) and the relation between changes in the bath correlation function and in the reduced dynamics and single-time averages of the system derived in Ref. (SpinBosonBounds, ), which give our approach a unique standing as a modeling technique based on an effective-environment concept but retaining the benefits of fully microscopic methods in terms of accuracy and rigor.

Our scheme maps a general open-system problem onto a Lindblad master equation for the system coupled to one or more small networks of interacting effective modes; the reduced dynamics of the system is then simulated by integrating this equation using standard numerical methods for Markovian problems and tracing out the oscillators. The surrogate modes are always at zero temperature regardless of the temperature of the original environment, giving the Lindblad equation a simple structure, and the interactions among them make a smaller number of modes necessary to account for the specific effect of any given environment than would be the case if they were all independent, with clear computational advantages. Not all modes need to be interacting; environmental spectral densities consisting of several terms may be reproduced using separate clusters of oscillators for each term, simplifying the calculation of their parameters while still exploiting the versatility of interacting oscillators in the rendering of each individual contribution.

As a first example of realistic use, we applied the method to a non-perturbative problem of a kind relevant for current research on coherent dynamics in biological molecular aggregates, obtaining accurate predictions (as confirmed by cross-checks with simulations performed using the well-established TEDOPA method) across the temperature range from absolute zero to room temperature with desktop-level computational resources. In addition, it was shown that by mapping non-Markovian problems to Markovian ones obeying a Lindblad quantum master equation, our approach can deal with short as well as long evolution times at comparable costs, making it a suitable tool for the simulation of long-lived dynamical features and relaxation to equilibrium. Next, we tested the performance of the technique on a much more challenging model inspired by the organic photovoltaics literature, demonstrating how the use of interacting surrogate modes allows for efficient simulation of systems so far only studied under much coarser approximations due to computational constraints. We have thus proved that the method is accurate, powerful and reliable and that there are classes of problems which no other approach known to us can successfully treat.

Future work on this project will be aimed mainly at improving the recipe for determining effective environment parameters, enhancing simulation efficiency and carrying out more case studies in order to better assess accuracy and performance, as well as applications to more systems of theoretical or experimental interest. Regarding the conceptual part of this work, we plan to undertake further analyses of the mapping problem from unitary to dissipative environments, in order to make the transformation more straightforward and possibly develop new TSO algorithms for different surrogate environment geometries. A deeper, system-dependent understanding of the error propagation from the environmental correlation function to the reduced dynamics could also be helpful in determining optimal figures of merit for individual problems, which would be a useful development for situations prioritizing accuracy over portability of the effective environment parameters. As a long-term goal, a fermionic extension of the entire approach based on the recent development of the relevant theoretical basis (Chen_FermionicTSO, ) is possible. Concerning the numerical implementations, it is our intention to refine and improve the codes for both the transformation and the simulations using various methods and possibly make them publicly accessible, as well as to investigate the possibility to integrate the two stages of our approach, for instance by choosing the surrogate environment geometry depending on the simulation routines for maximum efficiency.

Acknowledgements.

We wish to thank A. Mattioni, F. Caycedo-Soler, J. Lim, G. Gasbarri, R. Puebla and M. Paternostro for useful discussions, suggestions and feedback about the work presented in this paper. We further acknowledge support by the State of Baden-Württemberg through bwHPC for the use of the BwUniCluster and the German Research Foundation (DFG) through grant No. INST 40/467-1 FUGG for the use of the JUSTUS cluster. This work was supported by the ERC Synergy Grant BioQ, the EU Projects HYPERDIAMOND and AsteriQs, the BMBF projects DiaPol and NanoSpin, the Center for Integrated Quantum Science and Technology (IQST) and the Foundation BLANCEFLOR Boncompagni Ludovisi, née Bildt.

Appendix A Transformation to Surrogate Oscillators in detail

In Section III, we defined the general form of the effective environments used in our method, and sketched the transformation algorithm by which we obtain their parameters given a target correlation function. Here we will go through the procedure in detail, in order to give a clearer view of its more technical aspects.

A.1 Effective correlation function

The Hamiltonian for our effective oscillators in a chain configuration with hopping interactions is (with $\hbar=1$)

[TABLE]

and we consider a zero-temperature Lindblad dissipator

[TABLE]

acting locally on each mode. The interaction term with the system has the form

[TABLE]

where the interaction operators $F_{Rk}$ of the environment are linear in the creation and annihilation operators:

[TABLE]

Assuming factorizing initial conditions $\rho_{0}=\rho_{0S}\otimes\rho_{0R}$ with $\rho_{0R}=\bigotimes^{N}_{n=1}\ket{0}\!\bra{0}_{n}$, which is Gaussian and stationary under this dynamics, meaning it satisfies

[TABLE]

the correlation function

[TABLE]

is independent of the first evolution time $\tau$. We will drop the $\tau$ time argument from now on and also restrict our analysis to a single interaction operator ($m=1$ in Eq. (46)), so in the following the correlation function (48) will be denoted by $C^{R}(t)$. Writing it out explicitly in terms of the expression of $F_{R}$, we get

[TABLE]

since terms with two creation or annihilation operators and contributions proportional to $\langle b^{\dagger}_{m}(t)b_{n}(0)\rangle_{R}$ are zero for our initial vacuum state.

It is easy to show that the hopping coupling constants $g_{n}$ can be assumed real and positive without loss of generality in $C^{R}(t)$: define the canonical transformation

[TABLE]

for arbitrary real $\delta_{n}$. The creation operators $b^{\dagger}_{n}$ will transform with the opposite phase, preserving the canonical commutation relations. The free term in the Hamiltonian (44) and the dissipator (45) are invariant under this transformation; the hopping term in Eq. (44) and the interaction operator $F_{R}$ defined as in (47) are not:

[TABLE]

Taking $\delta_{n}$ such that $g_{n}e^{i(\delta_{n}-\delta_{n+1})}=|g_{n}|$, we may absorb the phase of the couplings in the still undetermined $c_{n}$, without restricting the physical picture in any way. Note that this leaves one of the $\delta_{n}$ still free as an overall phase in all the operator coefficients, which may be set e.g. so that $c_{1}$ or $c_{N}$ is real.

The free dynamics of the oscillators with no coupling to the system is given by the Lindblad equation

[TABLE]

Acting with the operator $b_{n}$ from the left on both sides and taking the trace, we get

[TABLE]

with

[TABLE]

where we have introduced the shorthand $\alpha_{n}\coloneqq-\frac{\Gamma_{n}}{2}-i\Omega_{n}$. Eq. (52) can be solved formally by diagonalizing the tridiagonal matrix $M$. Since $M$ is not Hermitian, one has

[TABLE]

with $\Lambda\coloneqq\mathrm{diag}(\lambda_{1},\dots,\lambda_{N})$ the diagonal matrix containing the eigenvalues, $S\coloneqq(\mathbf{u}^{1},\dots,\mathbf{u}^{N})$ a matrix made of arbitrarily normalized right eigenvectors $\mathbf{u}^{n}$ and $S^{-1}\coloneqq(\mathbf{v}^{1},\dots,\mathbf{v}^{N})^{T}$ its inverse, whose rows $(\mathbf{v}^{n})^{T}$ are left eigenvectors. Since $M$ is a symmetric matrix, left and right eigenvectors are the same, so $S^{-1}$ is just the transpose of $S$ up to normalization of the rows in such a way that

[TABLE]

Assuming that none of the eigenvalues are degenerate, which is always the case in numerical applications since the $\Lambda$ matrices with equal diagonal elements are a zero-measure set, the evolution of the expectation value $\langle b_{n}(t)\rangle_{R}$ is thus

[TABLE]

and extends to the correlation functions $\langle b_{n}(t)b^{\dagger}_{m}(0)\rangle_{R}$ by the quantum regression hypothesis, which is true by construction in the context of Theorem 1 [61]; since $\langle b_{n}(0)b^{\dagger}_{m}(0)\rangle_{R}=\delta_{nm}$ on our initial state, one has

[TABLE]

which can now be substituted into (49) to give the expression found in the main text:

[TABLE]

or

[TABLE]

in terms of the coefficients

[TABLE]

If degenerate eigenvalues $\lambda^{\mathrm{d}}_{k}$ are present, the time evolution in the corresponding subspace will be driven by $e^{\lambda^{\mathrm{d}}_{k}t}$ times growing powers of $t$; we did not consider this case for the sake of simplicity, but it may be useful to keep in mind that a mixed algebraic and exponential time dependence of correlation functions is not entirely ruled out by considering a Lindblad dynamics. If one wishes to explore this possibility in the TSO method, equality of two or more eigenvalues should be enforced at the level of the initial fit of the original correlation function $C^{E}(t)$ (see the next paragraph), since its spontaneous occurrence in the numerical procedure is virtually impossible.

A.2 Inversion problem from a given correlation function

To construct an effective environment whose $C^{R}(t)$ is as similar as possible to the $C^{E}(t)$ of a given unitary environment, we first fit $C^{E}(t)$ with a linear combination of $N$ complex exponentials $e^{\tilde{\lambda}_{n}t}$ weighted by complex coefficients $\tilde{w}_{n}$, with $N$ large enough to give an accurate fit, and then work backwards from Eq. (57) to find the parameters that give the best approximation of the target function.

Since the real parameters in $C^{R}(t)$ are $4N-1$ (taking into account the fact that $\sum^{N}_{n=1}w_{n}=\sum^{N}_{n=1}|c_{n}|^{2}$ is real and positive by construction) and it takes $5N-2$ real parameters ($N$ frequencies, $N$ damping rates, $N-1$ couplings and $N$ complex coefficients $c_{n}$ minus one overall redundant phase) to identify an effective environment, this is a highly nontrivial inversion problem, because the map from effective environments to correlation functions is both nonlinear and many-to-one. This means that existence or uniqueness of a solution to our problem are not guaranteed in general; furthermore, we must require $\Gamma_{n}>0$ for all $n$ in order for our effective master equation to be meaningful, which sets another important constraint.

It is useful to break down the problem into two parts: first an inverse eigenvalue problem leading from the $\tilde{\lambda}_{n}$ to the dynamical matrix $M$, and then a system of equations relating the coefficients $\tilde{w}_{n}$ to the interaction operator parameters $c_{n}$. This allows us to deal with the sign constraints on the rates once and for all in the first half of the solution procedure, and to exploit the fact that the $c_{n}$ only appear in the second.

To determine the relation between the eigenvalues and elements of the matrix $M$, it is not convenient to look for symbolic expressions for each eigenvalue in terms of the parameters, since these would necessarily involve high-degree roots of complex polynomials. A simpler approach is to consider the characteristic polynomial of $M$

[TABLE]

substitute the target eigenvalues $\tilde{\lambda}_{n}$ on the right-hand side and equate the coefficients of like powers of $\lambda$, which are geometrical invariants of any operator. The result is a system of equations of degrees 1 through $N$

[TABLE]

stating the invariance of the sums of principal minors order by order (the trace and the determinant appearing in the first and last equation being the simplest such invariants).

Now, Eq. (59) can be regarded as a parametric system of equations in the couplings $g_{n}$. With the $g_{n}$ fixed, it becomes an algebraic nonlinear system of $N$ equations in $N$ unknowns which can be solved numerically to give multiple sets of $\alpha_{n}$—i.e. frequencies $\Omega_{n}$ and rates $\Gamma_{n}$ whose sign can be checked directly—and therefore the entire dynamical matrix $M$.

Given a dynamical matrix $M$ obtained by choosing some set of $g_{n}$ and solving Eq. (59), its eigenvectors $\mathbf{u}^{n}$ and $\mathbf{v}^{n}$ can be substituted into the $w_{n}$ as defined in Eq. (58), which then become functions of the $c_{n}$ only and can be equated with the target values $\tilde{w}_{n}$

[TABLE]

to solve the second half of the problem. These $N$ complex equations are equivalent to $2N-1$ equations in $2N-1$ real unknowns, since the overall phase of all $c_{n}$ drops out of the left-hand side while on the right-hand side $\sum^{N}_{n=1}\tilde{w}_{n}=C^{E}(0)$ has no imaginary part.

A set of $c_{n}$ solving Eq. (60) does not always exist, so here we numerically minimize the Manhattan distance between the $w_{n}$ on the left-hand side and the $\tilde{w}_{n}$ instead. At this point, we have converted an arbitrary $(N-1)$-tuple of coupling constants into a trial effective correlation function

[TABLE]

which can be compared to the target $C^{E}(t)$ according to some figure of merit. We used the integral

[TABLE]

up to some final time $t_{\mathrm{max}}$ such that $C^{E}(t_{\mathrm{max}})\ll 1$ in all cases where we had a closed expression for it, and

[TABLE]

for some number of points $N_{\mathrm{max}}$ and timestep $\Delta t=t_{\mathrm{max}}/N_{\mathrm{max}}$ when $C^{E}(t)$ was only known in integral form and needed to be evaluated for each value of the time argument. This whole procedure can be carried out for many values of the couplings in the physical parameter region, ranking the corresponding trial correlation functions by their values of the figure of merit in search of an optimum, in the spirit of the error bounds in Ref. [47] which relate the absolute difference between correlation functions to the changes in the reduced dynamics.

To summarize the steps described above, in order to find an effective environment corresponding to some correlation function $C^{E}(t)$, we first fit it with complex exponentials, and then overcome the mismatch between the number of variables from this fit and the number of parameters in the effective environment by setting up a variational problem in the $g_{n}$ couplings between neighboring surrogate modes. We sample multiple $(N-1)$-tuples $(g_{1},\dots,g_{N-1})$ in a suitably sized open set $(0,g_{\mathrm{max}})^{N-1}$, solve Eq. (59) for each of them and then plug the eigenvectors of all physically acceptable matrices $M$ found into Eq. (60) to determine the $c_{n}$. The trial correlation functions $C^{R}_{\mathrm{trial}}(t)$ constructed from each set of parameters are ranked according to the estimators (61) or (62), depending on the original $C^{E}(t)$, and we search for the minimum of the figure of merit in the space of the $g_{n}$.

This variational problem is not convex in general: both the shape of the region in $g_{n}$-space leading to physically admissible solutions and the dependence of the cost functions defined in Eqs. (61) and (62) on the couplings can be highly nontrivial, with trenches, pointed features, local minima and gaps without any solutions at all appearing at unpredictable locations. We have also found no obvious patterns giving any hints as to the existence of a region of physically acceptable values of $g_{n}$ (for any $N$) such that there are exact solutions, or the form that such a region may have, based on the target parameters. All these mathematical features are very model-dependent; Fig. 12 shows some examples of parameter space shapes and cost function behavior for different correlation functions approximated using $N=3$ surrogate oscillators. Though none of the examples shown gave a $C^{R}(t)$ of sufficient accuracy for practical use due to the small number of effective modes, they nonetheless give a clear qualitative idea of the variety of possible outcomes. In general, sampling the parameter space efficiently is difficult, and we are looking for ways to improve this part of the algorithm.

Appendix B Special cases with exact solutions

We will now show some examples, both general and related to the specific systems treated in the main text, of analytical solutions of the inversion problem in specific cases.

B.1 One and two oscillators

The simplest possible effective environment is a single damped oscillator ($N=1$ in Eq. (44), with interaction operator $F_{R}=c(b+b^{\dagger})$ since the phase of $c$ can be set to zero). This yields a correlation function

[TABLE]

where the time dependence at $t<0$ is defined by $C^{R}(-t)\coloneqq C^{R*}(t)$ because the Lindblad equation only gives $C^{R}(t)$ for positive times, as discussed in the main text. The Fourier transform of this function is a Lorentzian of width $\Gamma$ centered in $\Omega$:

[TABLE]

A single sharp peak at zero temperature in the target spectral density can be mapped to a mode like this by simple nonlinear fitting of $C^{E}(t)$ with a complex exponential, as we did for the dimer simulations in the main text: in that case, the peaks were antisymmetrized Lorentzians so the frequency and damping rate of the effective mode matched those from the original spectral density almost exactly.

A less trivial, still exactly solvable case is given by two interacting oscillators and was already introduced by Garraway in Ref. [50] to show that not only sums but also differences of Lorentzians can be modeled by pseudomodes. In that paper, only one of the modes is coupled to the system (i.e. $c_{2}=0$); here we lift this assumption to show a more general result.

The general correlation function for $N>1$ has the form

[TABLE]

with $\Re[\lambda_{n}]<0$, where the fact that the $w_{n}$ are complex changes the form in the frequency domain from a simple linear combination of Lorentzians to

[TABLE]

thus adding a linear frequency dependence in the numerator of the term associated with each mode. In phenomenological approaches where there is no intention of accurately simulating a specific correlation function, one may use ad hoc combinations of weights and exponents to cancel terms in the numerator of the full $C^{R}(\omega)$ written as a single fractional polynomial and achieve a steeper fall-off in frequency than is possible with individual Lorentzians with positive coefficients (again, one example is given in Ref. [50]). Such strategies hardly generalize beyond specific applications but can be helpful to mitigate the error associated with the behavior of $C^{R}(\omega)$ near the origin, which does not comply with the fluctuation-dissipation theorem in general, as discussed in the main text.

In the case $N=2$, the eigenvalues $\lambda_{1,2}$ and weights $w_{1,2}$ depend on the effective environment parameters through the relations

[TABLE]

where $c_{2}$ is taken to be real by fixing the overall phase mentioned in the preceding section. The first equation is readily inverted parametrically in $g$:

[TABLE]

The domain of physically admissible solutions is the set of all $g$ such that $\Gamma_{1,2}=-2\Re[\alpha_{1,2}]>0$ and can be found by using the formula for the square root of a complex number $z=z_{\mathrm{R}}+iz_{\mathrm{I}}$

[TABLE]

Using Eq. (66) in the equation for the weights, this becomes parametric in $g$ as well:

[TABLE]

This equation, which relates the four real quantities $g$, $\Re[c_{1}]$, $\Im[c_{1}]$ and $c_{2}$ to the three real numbers determining $\tilde{w}_{1,2}$, may or may not have a solution depending on the $g$ chosen, as discussed earlier: if the solution exists only for $g$ outside the physical region in which $\Gamma_{1,2}=-2\Re[\alpha_{1,2}]>0$, then it is necessary to operate variationally and minimize the distance $D_{2}(g)\coloneqq\sum^{2}_{n=1}|\tilde{w}_{n}-w_{n}(g)|$.

If one assumes $c_{2}=0$ in (65), as is done in Ref. [50], then the whole system can be inverted explicitly:

[TABLE]

but the solution only exists if the $\tilde{\lambda}_{1,2}$, $\tilde{w}_{1,2}$ given are such that the expression on the right-hand side of the last equation has a vanishing imaginary part. This is because now the phase of $c_{1}$ has also decoupled from the problem, removing a second real degree of freedom and making the system overdetermined: the balance between equations and unknowns is thus restored by this real constraint appearing on the $\tilde{\lambda}_{1,2}$, $\tilde{w}_{1,2}$.

In our applications, we used pairs of effective modes to reproduce narrow antisymmetrized Lorentzians at nonzero temperature: we found no exact solution and had to minimize $D_{2}(g)$ in most cases, but e.g. for an antisymmetrized Lorentzian with peak frequency $\Omega_{AL}=215\,\mathrm{cm}^{-1}$ and width $\Gamma_{AL}=10\,\mathrm{cm}^{-1}$ at temperature $77\,\mathrm{K}$ the system can be solved exactly for $213.0\,\mathrm{cm}^{-1}<g<213.8\,\mathrm{cm}^{-1}$ (Fig. 13). Note that in all cases we considered, the best fit of such thermalized peaks with two modes was always obtained by mode frequencies close to zero and a strong coupling $g$ between the two; fitting the same function with two noninteracting modes at the positive and negative peak frequencies was consistently found to be a less accurate choice even for such seemingly obvious target functions.

One further possibility when dealing with very narrow high-frequency modes (so that the correlation function error around $\omega=0$ is very small) is to replace each such mode by a single oscillator, with $\Omega$ and $\Gamma$ equal to those of the antisymmetrized Lorentzian peak, and to initialize this mode in a Gibbs state at the bath temperature. Since independent oscillators or sets of oscillators have no correlations with each other in the initial state, this does not affect any other parts of the surrogate environment at hand in any way.

Such a mode would obey a full thermal Lindblad equation, with dissipator

[TABLE]

comprising emission and absorption terms with rates obeying detailed balance:

[TABLE]

where $n_{\Omega}(\beta)\coloneqq(e^{\beta\Omega}-1)^{-1}$ is the Bose-Einstein distribution. Coupling a thermalized mode to a system via a coefficient $c$ results in a correlation function $C^{R}(t)$ combining two exponential contributions with weights proportional to the emission and absorption coefficients, which translate to two Lorentzians centered at $\pm\Omega$ in the frequency domain:

[TABLE]

This transformation can be convenient when a single oscillator with a thermalized population requires a lower local dimension than a pair of coupled surrogate modes initialized in the vacuum would, since this would limit the memory requirements of the simulation. We used this method to account for the strongly coupled high-frequency mode of the polymer simulations in Section V of the paper.

B.2 Three oscillators

Adding a third oscillator, we found exact solutions for $c_{2}=0$, which we did not use in any of the simulations discussed in the main paper but can be useful in general.

For $N=3$, the system of eigenvalue equations is

[TABLE]

and one may remove $\alpha_{2}$ from the last two equations by using the first, so that $\alpha_{1}$ and $\alpha_{3}$ can be regarded as effective functions of the real parameters $g_{1}$ and $g_{2}$.

With $c_{2}$ set to zero, the whole inversion problem is determined, since the equations for the $\tilde{w}_{n}$ will determine the values of $g_{1}$ and $g_{2}$ instead. Setting the overall phase so that $c_{3}$ is real, the equations can be written as

[TABLE]

where the last line again features a real expression on the left-hand side and a complex one whose imaginary part must be zero on the right-hand side. Since the first equation is real by construction, there are five real equations in the five real variables $g_{1},g_{2},\Re[c_{1}],\Im[c_{1}],c_{3}$ in Eq. (73), so the existence of solutions is only subject to the constraint $\Gamma_{n}=-2\Re[\alpha_{n}]>0$.

If $c_{3}$ is also set to zero, then the system (73) becomes

[TABLE]

and can be inverted explicitly, giving $c^{2}_{1}$, $\alpha_{1}$ and $\alpha_{3}$. But now the system (72) is overdetermined: the trace gives $\alpha_{2}$, and the last two complex equations can give $g_{1,2}$ only if the $\tilde{\lambda}_{n}$ and $\tilde{w}_{n}$ happen to satisfy two real relations among themselves (one because $c_{3}$ was removed from the problem, another because the phase of $c_{1}$ is now irrelevant). In particular, the expressions whose imaginary part must vanish now appear on the right-hand side of the last two lines of the full solution

[TABLE]

Appendix C Effective parameters

We list here several sets of effective parameters used in the simulations discussed in the main text, along with the local dimensions of each mode in each set at convergence. The corresponding spectral densities are defined in Eq. (26), Eq. (33) and Eq. (34), respectively.