Efficient simulation of non-Markovian system-environment interaction

TEDOPA [38, 39] is a numerically exact and certifiable simulation method for open quantum system dynamics [46]. It applies to general systems under linear interaction with an environment modeled by a set of independent harmonic oscillators such that the total Hamiltonian H can be split into the system, the environment and the interaction between the two as

$$\begin{eqnarray}&&H={H}_{\mathrm{sys}}+{H}_{\mathrm{env}}+{H}_{\mathrm{int}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{env}}={\displaystyle \int }_{0}^{{x}_{\mathrm{max}}}{\rm{d}}x\;g(x){a}_{x}^{\dagger }{a}_{x},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{int}}={\displaystyle \int }_{0}^{{x}_{\mathrm{max}}}{\rm{d}}x\;h(x)({a}_{x}^{\dagger }+{a}_{x})A.\end{eqnarray} \tag{ 3 }$$

Here ${a}_{x}^{\dagger }$ and a_x denote the bosonic creation and annihilation operators corresponding to the environmental mode x. The coefficients $g(x)$ can be identified as the environmental dispersion relation. The interaction term ${H}_{\mathrm{int}}$ assigns each mode a coupling strength $h(x)$ between its displacement ${a}_{x}^{\dagger }+{a}_{x}$ and a general system operator A.

Together with the temperature, the functions $g(x)$ and $h(x)$ characterize the harmonic environment uniquely and define the spectral density$J(\omega )$ as

$$\begin{eqnarray}&&J(\omega )=\pi {h}^{2}[{g}^{-1}(\omega )]\displaystyle \frac{{\rm{d}}{g}^{-1}(\omega )}{{\rm{d}}\omega }.\end{eqnarray} \tag{ 4 }$$

Here ${g}^{-1}[g(x)]=g[{g}^{-1}(x)]=x$ and g(x) is monotonically growing. The quantity $\frac{{\rm{d}}{g}^{-1}(\omega )}{{\rm{d}}\omega }\delta \omega$ can be interpreted as the number of quantized modes with frequencies between ω and ω + δω (for δω → 0). We consider spectral densities subject to a hard cut-off at frequency ${\omega }_{\mathrm{hc}};$ this in turn defines the cut-off ${x}_{\mathrm{max}}$ in equations (2) and (3) as ${x}_{\mathrm{max}}={g}^{-1}({\omega }_{\mathrm{hc}})$.

TEDOPA uses a two-stage sequence to enable full treatment of the system and environment degrees of freedom. In a first step, an analytic transformation based on orthogonal polynomials converts the star-shaped system-environment structure into a one-dimensional geometry, where the system couples only to the first site of a semi-infinite chain of harmonic oscillators that contains only nearest-neighbour interactions. This is accomplished by use of a unitary transformation U, which defines new harmonic oscillators with creation and annihilation operators ${b}_{n}^{\dagger }$, b_n given by

$$\begin{eqnarray}&&{U}_{n}(x)=h(x){p}_{n}(x),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{b}_{n}^{\dagger }={\displaystyle \int }_{0}^{{x}_{\mathrm{max}}}{\rm{d}}x\;{U}_{n}(x){a}_{x}^{\dagger }.\end{eqnarray} \tag{ 6 }$$

Here ${p}_{n}(x)$ are orthogonal polynomials defined with respect to the measure ${\rm{d}}\mu (x)={h}^{2}(x){\rm{d}}x$ and the three-term recurrence relation

$$\begin{eqnarray}&&{p}_{k+1}(x)=(x-{\alpha }_{k}){p}_{k}(x)-{\beta }_{k}{p}_{k-1}(x),\end{eqnarray} \tag{ 7 }$$

with ${p}_{-1}(x)\equiv 0$ and k a positive integer or zero. This transformation can be performed analytically for specific spectral densities [39, 47]. For arbitrarily shaped spectral densities a numerically stable approach has been developed [38, 48]. Mappings using orthogonal polynomials have been shown to be exact for quadratic Hamiltonians [49].

The recurrence relation (7) results in the one-dimensional configuration with the Hamiltonian

$$\begin{eqnarray}\tilde{H} & = & {H}_{\mathrm{sys}}+{\eta }_{0}A({b}_{0}+{b}_{0}^{\dagger })+\displaystyle \sum _{n=0}^{\infty }{\omega }_{n}{b}_{n}^{\dagger }{b}_{n}\\ & & +\displaystyle \sum _{n=0}^{\infty }{\eta }_{n}({b}_{n}^{\dagger }{b}_{n+1}+{b}_{n}{b}_{n+1}^{\dagger }).\end{eqnarray} \tag{ 8 }$$

For a linear dispersion relation $g(x)={g}^{\prime }x$, ${\omega }_{n}={g}^{\prime }{\alpha }_{n}$ and ${\eta }_{n}={g}^{\prime }\sqrt{{\beta }_{n+1}}$. Due to the form of the emerging linear geometry, this first stage of TEDOPA is coined 'chain mapping'. One dimensional quantum many body systems can be efficiently simulated with the well established time dependent density matrix renormalization group (t-DMRG) algorithm [50–52]. Its application to evaluate the dynamics of the system and the transformed environment constitutes the second stage of TEDOPA. The long ranged correlations appearing in the original star-shaped geometry advise against the application of t-DMRG in that picture: it is the nearest neighbor structure that makes the numerical t-DMRG approach particularly efficient. Recent works consider the possibility to use generalized matrix product state formulations for treatment of star-shaped baths as well [53].

For a complete account of TEDOPA's inner workings refer to [38, 39, 54]. Suffice it to say that three main aspects distinguish it from other open quantum-system methods. First, it does not restrict the ratio $\lambda$ between inner-system couplings and system-environment couplings, unlike numerous other methods (which assume either $\lambda \gg 1$ or $\lambda \ll 1$). Second, the spectral density characterizing the system-environment interaction may assume any arbitrary shape, including a wide variety of sharp features that may be related to long-lived vibrational modes [38, 55, 56]. Such spectral densities are often encountered in spectral densities reconstructed from experimental data, e.g. in biological settings [57]. And third, while applicable to all temperatures, due to its scaling properties TEDOPA is inherently well-suited for simulations in the low-temperature domain.

Naturally, however, exact numerical methods tend to be costly and an effort to save on the associated computational demands is desirable. Where the transfer tensor method can be applied, these costs can be reduced and challenging long time simulations become accessible.

The transfer tensor method [41] reduces the numerical effort of TEDOPA simulations for a large class of non-Markovian environments. Its key idea is to relate the initial stages of the system's evolution to later times by efficiently determining the dynamical correlations built up between system and environment. This is achieved by the reconstruction of dynamical maps for short initial evolution times and their subsequent transformation into so-called transfer tensors.

A dynamical map is defined for an initial condition where the state of the system and the state of the environment are separable, and it fully determines the reduced state of the system $\rho (t)$ when applied to an initial condition $\rho (0)$

$$\begin{eqnarray}&&\rho (t)={ \mathcal E }(t,0)\rho (0).\end{eqnarray} \tag{ 9 }$$

For a time independent Hamiltonian H such as equation (1), equation (9) is related to the solution of the time-translationally invariant Nakajima–Zwanzig equation

$$\begin{eqnarray}&&\dot{\rho }(t)=-{\rm{i}}{{ \mathcal L }}_{s}\rho +{\displaystyle \int }_{0}^{t}{\rm{d}}{t}^{\prime }{ \mathcal K }(t-{t}^{\prime })\rho ({t}^{\prime }),\end{eqnarray} \tag{ 10 }$$

where ${{ \mathcal L }}_{s}\rho =[{H}_{\mathrm{sys}},\rho ]$ is the Liouvillian of the system alone and ${ \mathcal K }(t-{t}^{\prime })$, is the memory kernel arising from the system-bath interaction H_int.

The input for TTM is a set of dynamical maps ${{ \mathcal E }}_{k}={ \mathcal E }({t}_{k},0)$, where a discretisation time step $\delta t$ (${t}_{k}\equiv k\delta t$) is used. These dynamical maps are obtained easily from the evolution of all distinct initial basis states of the system's density matrix. The transfer tensors are then iteratively defined by

$$\begin{eqnarray}&&{T}_{n}={{ \mathcal E }}_{n}-\displaystyle \sum _{m=1}^{n-1}{T}_{n-m}{{ \mathcal E }}_{m}.\end{eqnarray} \tag{ 11 }$$

According to this definition, the transfer tensor T_k then quantifies the correlation in the dynamical map ${{ \mathcal E }}_{k}$ with the non-Markovian effects built up during the previous k time steps.

Further, the discretization of the memory kernel ${ \mathcal K }$ is directly related to these tensors by the time step $\delta t$

$$\begin{eqnarray}&&{T}_{k}={{ \mathcal K }}_{k}\delta {t}^{2},\end{eqnarray} \tag{ 12 }$$

where ${{ \mathcal K }}_{k}={ \mathcal K }({t}_{k})$ is the discretized memory kernel at time t_k. The system Hamiltonian is itself accessible from the first transfer tensor by

$$\begin{eqnarray}&&{T}_{1}=({\mathbb{1}}-{\rm{i}}{{ \mathcal L }}_{s}\delta t).\end{eqnarray} \tag{ 13 }$$

These tensors can be used to propagate the system to arbitrary later times as long as they cover the relevant part of the memory kernel. Assuming a finite coherence time of the bath ${t}_{\mathrm{bath}}$, the transfer tensors will decay sufficiently fast that a cutoff K can be defined such that T_k = 0 for $k\gt K$. Then, $\rho ({t}_{n})$ for $n\gt K$ can be expressed simply as

$$\begin{eqnarray}&&\rho ({t}_{n})=\displaystyle \sum _{k=1}^{K}{T}_{k}\rho ({t}_{n-k}).\end{eqnarray} \tag{ 14 }$$

TTM is applicable in a variety of interesting cases, scaling favorably in system size and length of the environment's correlation time. Due to its nature it also does not rely on assumptions about the system's parameters or environmental couplings. Thus it is usable as a powerful black box tool which, given initial trajectories, subsequently delivers the evolution trajectories for later times. For a more complete account of this tool refer to [41].

The combination of both methods facilitates the simulation of open quantum systems in regimes that were previously inaccessible. Exceptionally relevant is the ability to perform long-time simulations of low-temperature, highly structured harmonic environments at merely a fraction of the computational cost which would be necessary if only TEDOPA were applied. Evidence for this is provided by simple examination of some of the features of each of the methods. On the one hand, TEDOPA is based on a matrix product operator (MPO) description of the complete system plus environment density matrix. Settings leading to low occupation numbers of environmental oscillators—such as low temperatures—are especially suitable as they reduce the MPO's number and size. In addition TEDOPA is not limited to a specific analytical form of spectral density. On the other hand, recurrence effects originating at the chain's boundary limit the time for which accurate simulation is possible. Because TTM only requires sample system trajectories for as long as bath correlations are present, the required chain length is then not anymore determined by the target simulation time. The combination of TTM and TEDOPA is therefore most useful in highly non-Markovian regimes where bath correlation times are comparable to the maximum simulation time that TEDOPA can reach before recurrences appear.

Here we analyze relevant parametric cutoffs for the control of the accuracy of numerical simulations that combine TEDOPA and TTM.

In the case of TEDOPA one can distinguish between parameters related to the chain mapping and to the t-DMRG propagation. The semi-infinite chain of oscillators generated by the mapping necessarily requires the truncation of both the chain length N and the Hilbert space dimension d for each oscillator at a reasonable value. Those are native TEDOPA error sources and they have a direct effect on the recurrence time of the simulation and the maximum temperature that may be simulated, respectively. It should be noted that the error incurred by these two approximations can be upper-bounded rigorously by analytical expressions [46]. The effect of other relevant parameters, namely the MPO's matrix size $(\chi \times \chi )$ and the time step $\delta t$, are already well-known from the time-evolving block decimation (TEBD) algorithm [51]. To reach the accuracy required by TTM, care needs to be taken in adjusting these parameters to bound the total error of TEDOPA sufficiently. Some indications on how to accomplish this are provided in the present section.

The maximum time t_max before unphysical back-actions of the environment due to reflections at the end of the chain appear is related to the chain length N. Usually all chain coefficients are of the same order of magnitude and hence the simulation time  t_max scales roughly linearly with N. This reveals one of the benefits of the application of TTM on TEDOPA: the chain length can be truncated according to the length of the bath correlation time, allocating the simulation resources properly and shortening simulation times considerably in many cases of relevance (see figure 7). The exact relationship between N and t_max may be derived analytically through the use of Lieb–Robinson bounds [46] or numerically by trial-and-error: by setting the chain in an initial state $| 10\ldots 0\rangle$ and following the evolution of the number operator n on the first site, $O=n\otimes {\mathbb{1}}\otimes \ldots \otimes {\mathbb{1}}$, until a recurrence occurs.

The second native TEDOPA parameter is the local dimension d of the single oscillators constituting the environment. For a given temperature, the occupation of the single oscillators can be determined exactly, giving a rough scale of the necessary truncation level. Some error will be introduced necessarily but this can be upper-bounded analytically as explained in [46]. On the other hand, numerical benchmark calculations with increasing local dimensions will generally yield sufficiently accurate results.

A further subtlety in the chain mapping consists in the determination of an adequate hard cutoff frequency ${\omega }_{\mathrm{hc}}$ of the spectral density. For instance, the slow approach to zero for large frequencies of the Drude–Lorentz bath imposes a careful convergence check of the resulting physical behavior. For further discussions of these effects refer to [54].

While some error sources (like the cutoff in the chain length) introduce, if treated correctly, virtually no error at all, the matrix size $\chi$ necessarily does so due to the nature of the MPO. However, as already studied in the context of the TEBD algorithm, this error can be monitored during the time evolution [50]. This results in a quantity very similar to the discarded weight known from DMRG

$$\begin{eqnarray}&&{w}_{\mathrm{discarded}}=1-\displaystyle \sum _{i}{e}_{i}^{2}.\end{eqnarray} \tag{ 15 }$$

This quantifies the deviation from the targeted state using the discarded eigenvalues e_i. This error propagates in a non-trivial fashion and it is advisable to perform convergence checks in the dynamics under variation of the size of $\chi$. An additional source of error is derivated from the Suzuki–Trotter decomposition used in the TEBD part of TEDOPA.

It should be noted that the magnitude of the singular values kept during MPO-procedures should not fall below some threshold e₀. The transfer tensors determined by TTM do decay rapidly, falling to comparatively low magnitudes, and singular values corrupted by numerical noise deteriorates the interpretation of results as well as the propagation procedure.

Finally, for TTM the important quantity to keep track of in simulations is the norm of the memory kernel. This corresponds to the norm of the transfer tensors divided by the squared time step $\delta t$. This magnitude should exhibit a sufficiently fast decay so that the remaining tail can be neglected. Additionally, the time step $\delta t$ must be such that it provides a good resolution of the features of the memory kernel.

By construction, TEDOPA is inherently suited to treat spectral densities of arbitrary shape. When considering non-Ohmic spectral densities, Markovian approaches are well-known to anomalously suppress the effect of pure dephasing contributions [58, 59]. In this section we present an analysis of the dynamical effects of three different non-Ohmic spectral densities, namely

$$\begin{eqnarray}&&{J}_{1}(\omega )={\lambda }_{1}\;{\omega }^{3}\;{{\rm{e}}}^{-\omega /{\omega }_{c}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{J}_{2}(\omega )={\lambda }_{2}\;{\omega }^{5}\;{{\rm{e}}}^{-\omega /{\omega }_{c}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{J}_{3}(\omega )={\lambda }_{3}\;\sqrt{\omega }\;{{\rm{e}}}^{-\omega /{\omega }_{c}}.\end{eqnarray} \tag{ 20 }$$

To facilitate comparison, all of them exhibit the same exponential decay with ${\omega }_{c}=0.3\epsilon$ and are subject to a hard cutoff at ${\omega }_{\mathrm{hc}}\quad =\quad 10\epsilon$. Also the parameters ${\lambda }_{1}\quad =\quad 1.8\epsilon$, ${\lambda }_{2}\quad =\quad 1.0\epsilon$ and ${\lambda }_{3}\quad =\quad 0.6\epsilon$ have been chosen in such a way that they all share the same reorganization energy ${\lambda }_{r}=0.3\epsilon$. Thus the average interaction strength between system and environment is the same and the functional form of the spectral density is the factor responsible for disparate dynamics. The resulting dissimilar amplitudes and decay rates of the oscillations due to the different spectral densities are illustrated in figure 1. The tunneling strength is, as in previous section, ${\rm{\Delta }}\quad =\quad 0.6\epsilon$. One can observe that, for the fastest bath J₂, oscillations are sustained for a longer time, while this ability decreases for spectral densities centered in lower frequencies J₁ and almost disappears for the very slow bath represented by J₃. Some brief initial time is sufficient to generate the transfer tensors and predict the further evolution, whereupon high-accuracy TEDOPA simulations are used to verify these predictions.

Zoom In Zoom Out Reset image size

The suitability of the TEDOPA-TTM combination is supported by the fact that these simulations require on the order of just 100 tensors to converge to the exact results that have been obtained by full TEDOPA simulations, as shown in figure 2. This translates into about an order of magnitude faster results for TTM-TEDOPA combination than for TEDOPA alone. Further improvements in simulation speed are possible and are discussed in section 4.4.

Zoom In Zoom Out Reset image size

To further illustrate the power of our approach, we present results for a broad range of low to very low temperatures, up to $\beta \;=\;10\epsilon$. For the super-Ohmic spectral density ${J}_{2}(\omega )$ we show in figure 3 that it is possible to simulate the dynamics of a monomeric system at various inverse temperatures and the same system parameters as in the previous example. For the case of spectral density ${J}_{1}(\omega )$ we employ TTM to propagate the system until the steady state is reached (figure 4) and plot the steady-state occupation of the excited state for various inverse temperatures $\beta$ in figure 5. The insets in figures 3 and 4 show the memory kernel's decay over several orders of magnitude for the corresponding spectral densities. It is this decay which certifies the possibility to use the tensors for long-time propagation of the dynamics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The combination of TEDOPA and TTM is especially indicated for applications where accurate simulation of long time dynamics is crucial. The determination of absorption spectra belongs to this class of problems and we analyze here the more complex case of a dimeric system consisting of two coupled monomers.

The coupled dimeric system in the single excitation manifold consists of two excited states $| {e}_{1}\rangle$, $| {e}_{2}\rangle$ and a common ground state $| g\rangle$, and is described by the Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{sys}}={\epsilon }_{1}| {e}_{1}\rangle \langle {e}_{1}| +{\epsilon }_{2}| {e}_{2}\rangle \langle {e}_{2}| +J(| {e}_{1}\rangle \langle {e}_{2}| +| {e}_{2}\rangle \langle {e}_{1}| ),\end{eqnarray} \tag{ 21 }$$

where parameters ${\epsilon }_{1}\equiv \epsilon$, ${\epsilon }_{2}=2\epsilon$, and exchange interaction strength $J=0.6\epsilon$ are chosen. Each of the two systems is coupled to a bath. Without loss of generality we assume both environments are described by the same spectral density ${J}_{1}(\omega )$ and at temperature $\beta =\epsilon$.

The absorption spectrum is calculated as the Fourier transform of the two point correlation function of the dipole operator $\hat{\mu }={\mu }_{1}| {e}_{1}\rangle \langle g| +{\mu }_{2}| {e}_{2}\rangle \langle g| +{\rm{h.c.}}$

$$\begin{eqnarray}&&{C}_{\mu -\mu }(t,0)=\langle \hat{\mu }(t)\hat{\mu }(0)\rangle \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&=\;\mathrm{tr}[{{\rm{e}}}^{-{\rm{i}}{Ht}}\hat{\mu }{{\rm{e}}}^{{\rm{i}}{Ht}}\hat{\mu }\rho (0)],\end{eqnarray} \tag{ 23 }$$

between times t = 0 and $t=\tau$ such that the steady state has been reached at τ.

In the limit of weak interaction with the environment, the absorption spectrum emerging from Hamiltonian equation (21) exhibits two peaks in the one-exciton subspace. One of them is shown in the ${\lambda }_{1}=0.018\epsilon$ line (green) of figure 6, corresponding to the second excited state in the excitonic manifold at a wavelength of around $0.363\frac{c}{\epsilon }$. For higher coupling strengths with the environment, the emergence of the vibrational fine structure splits the peak in two, which is shown in the ${\lambda }_{1}=0.18\epsilon$ line (blue) and the ${\lambda }_{1}=1.8\epsilon$ line (black).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It will be interesting to compare the efficiency of the approach presented here with other approaches such as stochastic path integral methods which has recently been developed to calculate absorption and emission spectra [17] specifically for low temperatures and long times.

The ability of the TEDOPA-TTM combination to explore new simulation regimes is a consequence of the extraordinary savings in computational resources. We will explore these in terms of the 'wall time' t_w, the physical time required for the simulation to be executed as measured by an external clock.

Three factors have a direct influence on simulation time:

• bath coherence time ${t}_{\mathrm{bath}}$,
• chain length N and
• system dimension ${d}_{\mathrm{sys}}$.TTM requires the simulation of ${d}_{\mathrm{sys}}^{2}$ trajectories until ${t}_{\mathrm{bath}}$, one for each independent initial density matrix. Although this overhead may become inconvenient for systems of large dimension, the computation may be parallelized to avoid a scaling of t_w with ${d}_{\mathrm{sys}}$. Even without parallelization, numerical studies often require exploration of a large number of independent initial conditions anyway.

Due to the efficiency of multiplicative propagation with TTM (equation (14)), nearly the totality of the wall time t_w required for a simulation until ${t}_{\mathrm{sim}}$ is employed in the initial generation of the tensors until ${t}_{\mathrm{bath}}$ with TEDOPA. Therefore, one may consider t_w to be essentially independent of ${t}_{\mathrm{sim}}$. This makes the TTM-TEDOPA combination suitable for long time simulations, i.e. cases where ${t}_{\mathrm{sim}}\gg {t}_{\mathrm{bath}}$. There is an additional benefit in shortening simulations with TEDOPA to ${t}_{\mathrm{bath}}$, since this reduces the necessary chain length N.

The scaling of the wall time t_w necessary to perform a TEDOPA simulation of timestep $\delta t$ until ${t}_{\mathrm{bath}}$ can be expressed as

$$\begin{eqnarray}&&{t}_{w}\propto N\displaystyle \frac{{t}_{\mathrm{bath}}}{\delta t}\bar{t},\end{eqnarray} \tag{ 24 }$$

where dependence on three factors has been made explicit: the number of sites N, the number of time steps $\frac{{t}_{\mathrm{bath}}}{\delta t}$ and a factor $\bar{t}$ denoting the average wall time necessary to simulate one chain site during one time step $\delta t$. However, in order to avoid end-of-chain recurrences, for a simulation time ${t}_{\mathrm{bath}}$ one requires

$$\begin{eqnarray}&&N\propto {t}_{\mathrm{bath}}\cdot \bar{v},\end{eqnarray} \tag{ 25 }$$

sites N in the environment, given an average propagation speed $\bar{v}$ in the chain. Thus a total wall time of

$$\begin{eqnarray}&&{t}_{w}\propto \displaystyle \frac{\bar{v}\bar{t}}{\delta t}{t}_{\mathrm{bath}}^{2}\equiv c\cdot {t}_{\mathrm{bath}}^{2},\end{eqnarray} \tag{ 26 }$$

is needed where c is a scenario-dependent constant.

The global speedup provided by the TTM-TEDOPA combination is illustrated in figure 7 for three instances with different spectral densities. The near independence of t_w on ${t}_{\mathrm{sim}}$ is shown for large ${t}_{\mathrm{sim}}$. As shown in the inset, in reality t_w increases linearly with ${t}_{\mathrm{sim}}$, although with a negligible slope. For simulations with TEDOPA alone, the quadratic dependence expressed in equation (26) extends beyond ${t}_{\mathrm{bath}}$ until ${t}_{\mathrm{sim}}$.