Optimized auxiliary oscillators for the simulation of general open quantum systems

# Optimized auxiliary oscillators for the simulation of general open quantum systems

F. Mascherpa    A. Smirne Institut für Theoretische Physik, Universität Ulm, Albert-Einstein-Allee 11, D-89069 Ulm, Germany    D. Tamascelli Institut für Theoretische Physik, Universität Ulm, Albert-Einstein-Allee 11, D-89069 Ulm, Germany Università degli Studi di Milano, Dipartimento di Fisica, Via Celoria 16, I-20133 Milano, Italy    P. Fernandez Acebal    S. Donadi    S. F. Huelga    M. B. Plenio Institut für Theoretische Physik, Universität Ulm, Albert-Einstein-Allee 11, D-89069 Ulm, Germany
July 26, 2019
###### 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.

## 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, ).

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) (Diosi_QuantumStateDiffusion, ; Strunz_QuantumStateDiffusion, ; DiosiStrunz_QuantumStateDiffusion, ; Piilo_NonMarkovianQuantumJumpsPRA, ; Piilo_NonMarkovianQuantumJumpsPRL, ), 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 method (IlesSmith_ReactionCoordinate, ; IlesSmith_StructuredEnvironments, ) 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, ; 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 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 (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 an application of the method to a system in a structured environment relevant to current research investigating coherent effects in biomolecular aggregates; 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 (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:

 H\coloneqqHS⊗IE+IS⊗HE+HI (1)

with the free Hamiltonian of the environment

 HE\coloneqq∫∞0dωℏωa†ωaω

expressed in terms of creation and annihilation operators obeying the continuum canonical commutation relations , , and a general interaction term of the form (BreuerPetruccione, )

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

 ddtρ(t)=−iℏ[H,ρ(t)] (2)

and the initial state ; the reduced state of the system at time is obtained by taking the partial trace over the environmental degrees of freedom:

 ρS(t)=TrE[ρ(t)]. (3)

We are interested in the reduced dynamics of systems interacting with Gaussian environments, i.e. with linear in and and factorizing initial conditions with 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 and as given by the evolution of the oscillators with no coupling to the system:

 ⟨GEk(t)⟩E =TrE[U†E(t)GEkUE(t)ρ0E] (4) CEkk′(t+τ,τ) =TrE[U†E(t+τ)GEkUE(t)GEk′UE(τ)ρ0E] (5)

with .

### 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 Hamiltonian

 H′\coloneqqHS⊗IR+IS⊗HR+H′I (6)

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

 ddtρR(t)=LR[ρR(t)] (7)

where the new quantum Liouville superoperator

 LR[ρR]\coloneqq−iℏ[HR,ρR]+DR[ρR]

includes a dissipator

 DR[ρR]\coloneqqm∑i,j=1Λij(LRiρRL†Rj−12{L†RjLRi,ρR})

with a positive semidefinite rate matrix , which makes the dynamics non-unitary but ensures a completely positive and trace-preserving evolution at all positive times. The rate matrix , which we take to be constant in time, can always be brought into diagonal form by changing the basis of operators  (BreuerPetruccione, ); this is 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

 ddtρ(t)=L[ρ(t)] (8)

where

 L[ρ]\coloneqq−iℏ[H′,ρ]+D[ρ]

is the complete quantum Liouvillian for the system and the environment and embeds the dissipator into the full Hilbert space of the problem.

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

 HR\coloneqq∑nℏωnb†nbn

with and

 H′I\coloneqq∑lASl⊗FRl

with linear in the creation and annihilation operators, if one also takes the Lindblad operators linear in and and initial conditions with a Gaussian then the reduced dynamics of the system will only depend on the environment through and , again considering the free dynamics of the environment with no system attached, like in the unitary case:

 ⟨FRl(t)⟩R =TrR[FRleLRt[ρ0R]] (9) CRll′(t+τ,τ) =TrR[FRleLRt[FRl′eLRτ[ρ0R]]]. (10)

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 and correlation functions 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 (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

 ρUS(t)\coloneqqTrE[ρ(t)] (11)

for some system 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

 ρLS(t)\coloneqqTrR[ρ(t)] (12)

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 operators in and , with the corresponding and as well as the Lindblad operators of the non-unitary environment linear in the relevant creation and annihilation operators.

###### Theorem 1

(TSO_Theorem, ) Under the above assumptions, if

 ⟨FRk(t)⟩R=⟨GEk(t)⟩E∀k,t

and

 CRkk′(t+τ,τ)=CEkk′(t+τ,τ)∀k,k′,t,τ,

then

 ρLS(t)=ρUS(t)∀t.

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, since no equivalence theorem analogous to Theorem 1 for fermionic or other types of environments is known at the time of writing, we restrict our study to systems coupled to bosonic baths.

Finally, for physical reasons discussed in Section VI and thoroughly analyzed in (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, so we will apply the theorem in approximate form by looking for effective parameters such that and hence (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, and standard simulation methods for Lindblad equations may 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

 CEkk′(t+τ,τ)=CEkk′(t,0),

which will be denoted by in the following.

Any harmonic oscillator network obeying a Lindblad equation of the form (7), with quadratic in and , must also start from a stationary in order to give a time-homogeneous correlation function matrix for operators of the form

 FRk=∑n(cnkbn+c∗nkb†n). (13)

The condition for to be a stationary state of Eq. (7) is

 LR[ρ0R]=0. (14)

For the initial state of our effective environments, we will therefore need a Gaussian satisfying this property.

### iii.1 Ansatz and correlation function structure

The correlation functions of the auxiliary environment depend on all parameters appearing in , and the operators : 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 B.

We choose a free Hamiltonian corresponding to a chain of oscillators with a hopping interaction between nearest neighbors:

 HR\coloneqqN∑n=1ℏΩnb†nbn+N−1∑n=1ℏgn(bnb†n+1+b†nbn+1), (15)

where the couplings , as well as one of the coefficients in the interaction operators appearing in , can be assumed real without loss of generality if the are nonlocal, i.e. acting on all effective modes (see Appendix B). We complete the QME by adding local thermal dissipators at zero temperature acting on each oscillator:

 ddtρR(t)=−iℏ[HR,ρR(t)]+N∑n=1Γn(bnρR(t)b†n−12{b†nbn,ρR(t)}) (16)

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

 ρ0R=N⨂n=1|0⟩⟨0|n. (17)

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 and related by Hermitian conjugation. For one has

 ddt⟨bn(t)⟩R=N∑m=1Mnm⟨bm(t)⟩R, (18)

where

 Mnm:=(−Γn2−iΩn)δnm−i(gmδnm+1+gm−1δnm−1). (19)

The two-time correlation function 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 obey the same equations of motion as the single-time expectation values  (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 and plugging the result (as well as its conjugate , which is identically zero for our initial state (17)) into the expression of in terms of the operators as given in Eq. (13), one finds

 CRkk′(t)=N∑n=1(wn)kk′eλnt (20)

where are the eigenvalues of the matrix defined in Eq. (19), which we assume to be non-degenerate for simplicity (see Appendix B for further discussions), and

 (wn)kk′=N∑l,m=1clkc∗mk′unlvnm (21)

are complex coefficients obtained from the definition of the operators and the left and right eigenvectors and corresponding to each , normalized in such a way that as discussed in Appendix B. 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 . Our goal is to find the matrix elements and operator coefficients of some operator as given in Eq. (13) such that the resulting effective correlation function is as close as possible to .

The form of in terms of and is given by Eq. (20), where the eigenvalues and weights can be thought of as functions of the free parameters , , and—only the with , where is the number of oscillators making up the effective bath.

In order to determine the values of these free parameters such that , we proceed in two steps. First, we perform a nonlinear fit on using damped exponentials with complex coefficients

 CE(t)⟶~CE(t)=N∑n=1~wne~λnt, (22)

for instance using Prony analysis (Marple_Prony, ).

Then we solve the problem of matching or getting as close as possible to the and with the and 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 ( 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 need to be taken into account as well, so it is typically necessary to minimize some distance between and instead of exactly matching the best fit ; 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 , it is possible to invert the equations exactly and obtain valid effective bath parameters; we have listed a few explicit solutions in Appendix C and will put some of them to use in our example applications. For general , we devised a variational recipe to carry out our Transformation to Surrogate Oscillators (TSO) systematically. This is also described in detail in Appendix B and can be summarized as follows:

• Sample random points in a suitably sized open set , to be used as coupling constants.

• Substitute each -tuple into the equations relating the complex eigenvalues to the : this will give rates and frequencies such that the eigenvalues match; keep only the solutions with all .

• Compute the left and right eigenvectors of the matrices corresponding to each solution found, plug them into Eq. (21) and minimize a distance (e.g. the Manhattan distance ) between and by varying the .

• Assess overall accuracy of the solutions found and rank the corresponding according to a meaningful figure of merit, such as the integral from (SpinBosonBounds, ).

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 : 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 . 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.

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 —or equivalently, the underlying spectral density of the unitary environment—as shown in Fig. 3. Note that decoupling all oscillators from each other, i.e. taking all (which corresponds to requiring all to be real and positive in Eq. (22)), one recovers the noninteracting pseudomodes of (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 through the position operator of each oscillator (for a more succinct notation, we will leave the tensor products implicit and use natural units , from now on):

 H=HS+∫∞0dωωa†ωaω+AS∫∞0dωg(ω)(aω+a†ω). (23)

This type of coupling for microscopic models is one of the most common in the OQS literature (BreuerPetruccione, ; Weiss, ; CaldeiraLeggettModel, ; Leggett_SpinBoson, ; DeVegaAlonso_NonMarkovianity, ). The correlation function of the interaction operator for a thermal initial state at inverse temperature is

 (24)

where the spectral density is related to the frequency-dependent coupling strength through

 J(ω)=πg2(ω)

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 best approximating their behavior near the origin, where they are always zero; a is called Ohmic if , and super-(sub-)Ohmic if (). The spectral density and temperature uniquely determine 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

 CEβ(ω)=∫∞−∞dtCEβ(t)eiωt=(1+coth(βω2))(J(ω)θ(ω)−J(−ω)θ(−ω)), (25)

where is the Heaviside step function, is always real; at temperature , it is just . In fact, may itself be regarded as a spectral density defined over a new environment, which encompasses both positive- and negative-frequency modes and gives the correlation function if initialized in the vacuum (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 are only defined at positive times; we extend the definition to negative times by imposing the same symmetry

 CR(−t)\coloneqqCR∗(t)

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 with an exponential cutoff:

 J(ω)=πωe−ωΩc. (26)

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 corresponding to the spectral density (26) can be determined analytically as

 CEβ(t)=Ω2c(1+iΩct)2+1β2(ψ′(1+1+iΩctβΩc)+ψ′(1+1−iΩctβΩc)) (27)

where is the polygamma function of order one.

Performing our TSO on this correlation function at temperature according to the recipe described in the previous subsection, for we determined the parameters given in Table 1; Fig. 3 shows the Fourier-transformed effective correlation function obtained using these parameters and the target for comparison. As can be seen from the plot, four interacting oscillators were enough to obtain a very accurate , with a peak in the error around 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 is hardest to match: this is because any is analytical around zero by construction, whereas 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

 H=ω02σz+∫∞0dωωa†ωaω+k2σz∫∞0dωg(ω)(aω+a†ω) (28)

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 and 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 (BreuerPetruccione, ) for a derivation)

 ρ01(t)=ρ∗10(t)=e−iω0t+k2Γ(t)ρ01(0) (29)

with

 Γ(t)\coloneqq∫∞0dωπJ(ω)coth(βω2)cos(ωt)−1ω2.

Using the cutoff frequency of the environment as our energy scale, we set the parameter values and , 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 , with in terms of the eigenstates of , and we simulated the reduced dynamics at three different temperatures , and . 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 for the first two cases, and with for the high-temperature regime; the parameters are given in Appendix A, 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 at all times and temperatures (the populations and are both equal to throughout the evolution), as the overlap between the numerical (solid lines) and exact (dashed lines) solutions shows. The pure quantum decoherence at induces an algebraic decay asymptotically proportional to , while at the damping becomes exponential; a stronger effective coupling regime, which is determined both by and the strength of thermal effects, induces faster relaxation in the system dynamics.

The plots in the insets show the error figure

 Ef(t)\coloneqq|f(t)−fNum(t)||f(t)|+|fNum(t)|, (30)

which is identical for (solid line) and (dashed line). This is a better estimator for the accuracy of the simulations than e.g. the absolute difference 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 , 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 levels for most oscillators and a maximum of 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 A physically relevant application

Few-body systems non-perturbatively coupled to environments with structured spectra are ubiquitous in many fields ranging from biological physics and chemistry (HuelgaPlenio_QBio, ; Pelzer_Transport, ) to condensed matter (RibeiroVieira_Transport, ; Haase_Metrology, ) and 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 models and numerical methods for non-Markovian OQS. As a first example of an application of our simulation method to a system relevant for current research, we will now show some results for experimentally measurable optical properties in a model inspired by research in coherent charge and energy transfer in biological molecular aggregates (Scholes_QBioNature, ).

### v.1 The model

We consider 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 simulated the system and computed its absorption spectrum at temperatures ranging from to , comparing the results for two different versions of the environmental spectral density in order to identify the optical signatures setting them apart: in particular, we sought to determine the differences between absorption spectra obtained in the presence or absence of a strongly coupled vibrational mode in addition to a broad background spectral density.

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

 (31)

where the two monomers have on-site energies and and interact through a hopping coefficient , and is the state with both monomers excited. We assume our system to be probed by a low-intensity laser source, as it is in most absorption experiments (MayKuehn, ), and will therefore focus only on optical transitions between the ground state and the single-excitation manifold spanned by the states , ignoring the doubly excited state in the following since its contribution is typically negligible. Then, setting as our reference energy, we are left with an effective two-level Hamiltonian for the single-excitation manifold

 H1ex=2∑n=1En|En⟩⟨En|+J(|E1⟩⟨E2|+|E2⟩⟨E1|), (32)

whose eigenstates have an energy gap of .

The local excited states 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 the left panel of Fig. 5.

We first studied the problem for a spectral density in the super-Ohmic form introduced by Adolphs and Renger (AdolphsRenger, )

 JAR(ω)\coloneqqπ2⋅9!2∑a=1ρaω5Ω4ARae−√ω/ΩARa, (33)

where the two cutoff frequencies are and the weights of the two terms are . This spectral density accounts for environmental noise across a broad frequency range and will be referred to as in the following. Then we repeated the analysis for a spectral density featuring a strongly coupled vibrational mode in addition to the broad background: the mode was represented by adding an antisymmetrized Lorentzian peak

 JAL(Ω,Γ,S;ω)\coloneqqS8ΓΩ(4Ω2+Γ2)ω(4(ω−Ω)2+Γ2)(4(ω+Ω)2+Γ2), (34)

to the Adolphs–Renger background, and we chose a frequency , near resonance with the system, a width corresponding to a decay time , and a Huang–Rhys factor placing it in a strong-coupling regime with the system. We will refer to the spectral density with the strongly coupled mode as . The reorganization energies corresponding to the two environments are

 λ0 =∫∞0dωπJ0(ω)ω=2∑a=1ρaΩARa=19.93cm−1 λ1 =∫∞0dωπJ1(ω)ω=λ0+SΩ=69.93cm−1

and a plot of both spectral densities is given in Fig. 6.

The total Hamiltonian of our problem

 Htot=H1ex+2∑n=1∫∞0dωn(ωna†ωnaωn+g(ωn)|En⟩⟨En|(aωn+a†ωn)), (35)

can be rewritten in terms of the ‘common-mode’ and ‘relative’ creation and annihilation operators parametrized by a single frequency and . The common-mode terms then decouple from the system and the Hamiltonian reduces to

 H=H1ex+∫∞0dωωa†ωaω+1√2(|E1⟩⟨E1|−|E2⟩⟨E2|)∫∞0dωg(ω)(aω+a†ω) (36)

in terms of the ‘relative’ modes only. Note that the ground state is decoupled from the single-excitation subspace, since no off-diagonal terms involving it appear in the Hamiltonian (36). The dynamics thus 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 the right panel of Fig. 5.

Absorption experiments probe the linear response of the system; the light from the laser source can be described as interacting with the local dipole moment operators , where is the classical dipole moment of the -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 over the initial stationary state

 ρ0Abs\coloneqq|g⟩⟨g|ρβ, (37)

where the bath is in a thermal state at inverse temperature 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 , applied at times and : in terms of the overall unitary evolution, one has

 Cμ(t)\coloneqqTr[U†(t)(→μ†1+→μ†2)U(t)⋅(→μ1+→μ2)ρ0Abs]. (38)

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 , 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.

We set the ansatz for the geometry of the dimer in order to simplify the form of the dipole correlation function. Normalizing its value to , becomes

 (39)

The absorption spectrum is then given by

 SAbs(ω)\coloneqqωIlimtmax→∞∫tmax0dtiCμ(t)eiωt. (40)

### v.2 TSO, results and simulation metrics

In order to calculate and hence for our model, we determined effective parameters for the two spectral densities and at temperatures , () and (), performing the TSO separately on the Adolphs–Renger background (33) and the antisymmetrized Lorentzian peak (34). This corresponds to assigning a separate effective environment to each additive part of the spectral density and can be a convenient strategy to break down structured spectra, as mentioned in an earlier section and shown in Fig. 3. The effective parameters for the two environments at each temperature are then given by a common set accounting for the background and the extra parameters reproducing the peak; the Adolphs–Renger correlation function required oscillators at all three temperatures, and the Lorentzian mode was replaced by one effective oscillator at and two interacting ones at using the exact methods for described in Appendix C. All parameters of the effective environments are given in Appendix A.

The dipole correlation function can be calculated by solving the effective Lindblad equation with the initial pseudo-state . Since the memory required for a direct integration of the 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 for sparse Lindblad superoperators such as ours, while a master equation integrator requires at least ). For correlation functions, a convenient decomposition taking our initial condition into account is described in (DalibardCastinMolmer_MCWF, ): it can be shown that