Efficient Real-Time Path Integrals for Non-Markovian Spin-Boson Models

TL;DR

This paper introduces a new numerical method for simulating non-Markovian spin-boson models using real-time path integrals, accurately capturing dynamics including revivals caused by environmental memory effects.

Contribution

The authors reformulate finite memory approximations in path integral simulations to exactly reproduce non-Markovian dynamics, enabling more accurate modeling of complex quantum systems.

Findings

New method accurately reproduces exact dynamics

Reveals population revivals due to non-Markovian effects

Enables simulation of previously intractable parameter regimes

Abstract

Strong coupling between a system and its environment leads to the emergence of non-Markovian dynamics, which cannot be described by a time-local master equation. One way to capture such dynamics is to use numerical real-time path integrals, where assuming a finite bath memory time enables manageable simulation scaling. However, by comparing to the exactly soluble independent boson model, we show that the presence of transient negative decay rates in the exact dynamics can result in simulations with unphysical exponential growth of density matrix elements when the finite memory approximation is used. We therefore reformulate this approximation in such a way that the exact dynamics are reproduced identically and then apply our new method to the spin-boson model with superohmic environmental coupling, commonly used to model phonon environments, but which cannot be solved exactly. Our new method allows us to easily access parameter regimes where we find revivals in population dynamics which are due to non-Markovian backflow of information from the bath to the system.

Equations

$$H=$$

$$=$$

$${\hat{\rho }}_R$$

$$\equiv \sum _S{\rho }_S|S⟩⟩\equiv |{\rho }_R⟩⟩,$$

$$|{\rho }_R(t)⟩⟩={\mathrm{Tr}}_B\left[{\text{e}}^{\mathcal{L}t}|{\rho }_R(0)⟩⟩|{\rho }_B⟩⟩\right],$$

$${\text{e}}^{\mathcal{L}\Delta t}\approx {\text{e}}^{{\mathcal{L}}_B\Delta t}{\text{e}}^{{\mathcal{L}}_0\Delta t},$$

$$\begin{gathered} ⟨⟨S_N⟩{\rho }_R(t_N)⟩=\sum _{S_0\dots S_{N-1}}F\left(\{S_k\}\right)I\left(\{S_k\}\right) \\ \times ⟨⟨S_0⟩{\rho }_R(0)⟩. \end{gathered}$$

$$F\left(\{S_k\}\right)=$$

$$=$$

$$I\left(\{S_k\}\right)=$$

$$\eta_{k-k^{\prime}}=\begin{dcases*}\int_{t_{k-1}}^{t_{k}}\int_{t_{k^{\prime}-1}}^{t_{k^{\prime}}}C(t^{\prime}-t^{\prime\prime})dt^{\prime\prime}dt^{\prime}&\text{ $k\neq k^{\prime}$}\\ \int_{t_{k-1}}^{t_{k}}\int_{t_{k-1}}^{t^{\prime}}C(t^{\prime}-t^{\prime\prime})dt^{\prime\prime}dt^{\prime}&\text{ $k=k^{\prime}$}\end{dcases*},$$

$$C(t)$$

$$={\int }_0^{\infty }d\omega J(\omega )(\coth (\omega /2T)\cos (\omega t)-i\sin (\omega t)),$$

$$|{\rho }_R^A(t+{\tau }_c)⟩⟩={\text{e}}^{{\mathcal{L}}^A{\tau }_c}|{\rho }_R^A(t)⟩⟩,$$

$$|S_k^A⟩⟩=|S_k⟩⟩|S_{k-1}⟩⟩\dots |S_{k-\Delta k+1}⟩⟩,$$

$$|{\rho }_R(t_k)⟩⟩=\sum _{S_{k-1}\dots S_{k-\Delta k+1}}⟨⟨S_{k-1}|\dots ⟨⟨S_{k-\Delta k+1}⟩{\rho }_R^A(t_k)⟩.$$

$$⟨⟨S_k|{\text{e}}^{{\mathcal{L}}_0\Delta t}|S_{k-1}⟩⟩=⟨⟨S_k⟩S_{k-1}⟩⟩={\delta }_{s_k^{+}{s}_{k-1}^{+}}{\delta }_{s_k^{-}{s}_{k-1}^{-}}.$$

$$⟨⟨S_N⟩{\rho }_R(t_N)⟩={\text{e}}^{\Gamma (S_N,t_N)}⟨⟨S_N⟩{\rho }_0⟩,$$

$$\begin{gathered} \Gamma (S_N,t_N)=-(s_N^{+}-s_N^{-}{)}^2\mathrm{Re}[\eta (t_N)] \\ -i\left({\left(s_N^{+}\right)}^2-{\left(s_N^{-}\right)}^2\right)\mathrm{Im}[\eta (t_N)]. \end{gathered}$$

$$\eta (t_N)=\sum _{k=0}^N\sum _{k^{\prime }=0}^k{\eta }_{k-k^{\prime }}={\int }_0^{t_N}{\int }_0^{t^{\prime }}C(t^{\prime }-t^{\prime \prime })d{t}^{\prime \prime }d{t}^{\prime },$$

$$\eta_{k-k^{\prime}}=\begin{dcases*}\eta(t_{k-k^{\prime}+1})-2\eta(t_{k-k^{\prime}})+\eta(t_{k-k^{\prime}-1})&\text{ $k\neq k^{\prime}$}\\ \eta(t_{1})&\text{ $k=k^{\prime}$}\end{dcases*}.$$

$$\eta (t)\approx \eta ({\tau }_c)+(t-{\tau }_c)\frac{\eta ({\tau }_c+\Delta t)-\eta ({\tau }_c)}{\Delta t}.$$

$$\eta (t)=$$

$$\approx$$

$$=$$

$$J(\omega )=\frac{\alpha }{2}\frac{{\omega }^{\nu }}{{\omega }_c^{\nu -1}}{\text{e}}^{-\omega /{\omega }_c},$$

$$\sum _{k=0}^N\sum _{k^{\prime }=k-\Delta k}^k{\tilde{\eta }}_{k-k^{\prime }}=\sum _{k=0}^N\sum _{k^{\prime }=0}^k{\eta }_{k-k^{\prime }}=\eta (t),$$

$$\tilde{\eta}_{k-k^{\prime}}=\begin{dcases*}\eta_{k-k^{\prime}}+\sum_{j=0}^{k-\Delta k-1}\eta_{k-j}\equiv\eta_{\Delta k}^{k}&\text{ $k-k^{\prime}=\Delta k$}\\ \eta_{k-k^{\prime}}&\text{ $k-k^{\prime}<\Delta k$}\end{dcases*}.$$

$${\eta }_{\Delta k}^k$$

$$=\eta (t_k)-\eta (t_{k-1})-\eta ({\tau }_c)+\eta ({\tau }_c-\Delta t).$$

$$D({\rho }_1(t),{\rho }_2(t))=\frac{1}{2}\mathrm{Tr}|{\rho }_1(t)-{\rho }_2(t)|,$$

$$\Delta I(t)=\frac{d}{dt}D({\rho }_1(t),{\rho }_2(t)),$$

$${\alpha }_z(t)=A{e}^{-{\gamma }_{1}t}\cos ({\omega }_{1}t)+B{e}^{-{\gamma }_{2}t}\sin ({\omega }_{2}t+\varphi ),$$

$$\begin{gathered} F\left(\{S_k\}\right)I\left(\{S_k\}\right)=\left(\prod _{j=1}^N\prod _{j^{\prime }=1}^j\tilde{I}(S_j,S_{j^{\prime }})\right) \\ \times ⟨⟨S_1|e^{{\mathcal{L}}_0\Delta t}|S_0⟩⟩ \end{gathered}$$

Full text

I Introduction

Finding accurate descriptions of open quantum systems strongly coupled to external environments is essential for understanding how quantum systems lose their coherence. Chirolli and Burkard (2008); Breuer and Petruccione (2002). For example, the behaviour of quantum dots in semiconductors Madsen et al. (2011); Ramsay et al. (2010a, b); Wilson-Rae and Imamoğlu (2002); McCutcheon (2016) and NV centres in diamonds Stanwix et al. (2010); Wang and Takahashi (2013) can both require a strong coupling description. Strong coupling inevitably leads to the emergence of non-Markovian phenomena in such systems, and this has been experimentally demonstrated Liu et al. (2011), opening up the potential for exploiting non-Markovianity as a resource in developing quantum technologies Bylicka et al. (2014); Rivas et al. (2010).

In this paper we will introduce a new way of modelling a strongly coupled environment, developing previous approaches based on Feynman’s path integral formulation of open quantum systems Feynman and Vernon (1963); Feynman and Hibbs (1965); Caldeira and Leggett (1983); Weiss (2012). In particular, we will show how a technique based on discretisation of the Feynman influence functional – the so-called Augmented Density Tensor (ADT) can be modified to significantly improve the convergence of simulations in its numerical implementation. This allows us to study the spin-boson model in a very strong coupling regime that shows clear non-Markovian behaviour that we quantify with the widely-used trace distance measure Breuer et al. (2009).

When the system-environment coupling is weak, Born-Markov master equations Lindblad (1976); Breuer and Petruccione (2002) provide a perturbative approach in which the environment is assumed to be Markovian i.e. memoryless. While this approach is very successful within its range of applicability Breuer and Petruccione (2002); Nazir and Mccutcheon (2016), for many physical applications the (rather severe) approximations made when deriving these types of master equation are not justified. A common approach to go beyond the weak coupling regime is to use a polaron-transformed master equation McCutcheon and Nazir (2011); Roy and Hughes (2012); Pollock et al. (2013); Nazir and Mccutcheon (2016); Jang et al. (2008); Jang (2009). Such an approach has been used to understand features which could not be explained in weak coupling for semiconducting quantum dots McCutcheon and Nazir (2011), circuit QED Marthaler et al. (2011), energy transfer in biological systems Kolli et al. (2012); Zimanyi and Silbey (2012); Pollock et al. (2013) and Bose-Einstein condensation of photons in optical microcavities Kirton and Keeling (2013, 2015). However, a polaron master equation comes at the expense of introducing a restriction on the renormalized system Hamiltonian terms. Other master equation based techniques such as Nakajima-Zwanzig projection operator equations Breuer and Petruccione (2002), time-convolutionless master equations Breuer and Petruccione (2002), reaction co-ordinate methods Iles-Smith et al. (2014) and Keldysh-Lindbald equations Müller and Stace (2017) are able to pertubatively go beyond both the Born and Markov approximations but are still, in practice, limited to restricted parameter regimes.

One can overcome this restriction by using non-perturbative approaches. One class of such methods uses the Feynman path integral representation and formally integrates out the environmental degrees of freedom. Then all effects of the system-environment coupling are described by an influence functional Feynman and Vernon (1963) that acts only on the reduced system trajectories. This representation of an open quantum system has proved useful in developing both analytical and numerical methods. As an example, the non-interacting blip approximation Caldeira and Leggett (1983) is a successful non-perturbative analytical method for understanding the spin-boson model, derived from the influence functional. However, this is only applicable for a relatively small range of parameters.

A versatile numerical method for finding the dynamics of these kinds of systems is the ADT scheme, first introduced as the Quasi-Adiabatic Propagator Path Integral (QUAPI) by Makri and Makarov (1995a, b). This is, in principle, applicable to environments with arbitrary spectral density. The main drawback to this technique is the exponential scaling of storage requirements with system size, though in some cases ways of reducing these requirements are known Sim (2001); Makri (2017). Hence, the primary use of this method is to calculate the of dynamics of few-level systems in contact with a bosonic reservoir. The ADT has been used to calculate equilibrium correlation functions Shao and Makri (2002) and system steady states Makri and Makarov (1995a). In addition, multiple spatially separated sites coupled to the same bath can be accounted for Nalbach et al. (2010), and additional Markovian dynamics in the reduced system can be incorporated at no additional computational cost Barth et al. (2016). Practically, the ADT method has been put to use effectively in, for example, benchmarking master equation methods Chang et al. (2013); McCutcheon et al. (2011), simulating systems in the difficult-to-access regime of subohmic system-bath coupling Nalbach and Thorwart (2010) and investigating dissipative dynamics of charge qubits in realistic environments Thorwart et al. (2005); Eckel et al. (2006). Very recently, it was shown that, as long as the operator which couples to the environment only acts in a small section of Hilbert space much larger systems can be treated at only small numerical cost Cygorek et al. (2017).

Other numerical approaches are able to accurately describe non-perturbative dynamics in spin-boson systems. The hierarchical equations of motion method Tanimura and Kubo (1989) is regularly used to benchmark master equation approaches Iles-Smith et al. (2014). Techniques based on the numerical renormalisation group have been successfully applied to quantum impurity dynamics by mapping baths characterized by particular spectral density forms into semi-infinite one-dimensional chains Anders and Schiller (2005); Bulla et al. (2005); Anders et al. (2007); Bulla et al. (2008). Ansatz wavefunctions such as matrix product states Yao et al. (2013) can also be used following a similar mapping Prior et al. (2010), while approaches based on a Monte Carlo sampling of the path integral can be used to look at dynamics Egger and Mak (1994); Kast and Ankerhold (2013). For a recent review of many of the techniques mentioned above see de Vega and Alonso (2017).

The key approximation made in the ADT scheme is that the bath is assumed to have a sharply defined finite memory time, and so non-Markovian effects are given finite range. This is intuitively justified since, for infinite bosonic baths, all correlations decay to negligible size in finite time. The performance of this method has been successfully tested against other numerical methods Nalbach et al. (2011) and exact analytics Thorwart et al. (2000). However, the full consequences of the sharp memory time cutoff and how this affects convergence have not been addressed. Moreover, it has been recognized that for superohmic environments using this approximation results in both quantitatively and qualitatively incorrect asymptotic long time behaviour Vagov et al. (2011); G Nahri1 et al. (2016). In this paper we will show that, in certain cases, throwing away small long time correlations beyond the cutoff can have dramatic effects, including unbounded growth of density matrix elements. We then propose a less severe way to make this approximation which does not suffer from the same problems.

The structure of this paper is as follows: In Sec. II we introduce the ADT scheme in detail. Following this, in Sec. III, we investigate the results of enforcing the finite memory approximation upon the exact solution of a pure dephasing model and identify the class of bath spectral densities for which the finite memory approximation has most drastic effects. We then propose an alternative way of taking the finite memory approximation such that the exact solution of the dephasing model is reproduced at all times for arbitrary spectral densities. In Sec. IV we go on to apply our method to the spin-boson model with a highly non-Markovian superohmic bath. We find that our modified memory cutoff significantly improves convergence in this model and allows us to observe revivals in the population dynamics due to excitation exchange with the environment. We also calculate the trace distance measure of non-Markovianity in this model, finding regimes of non-Markovian behaviour that we compare to the population revivals. Finally, in Sec. V we summarise our results.

II ADT Scheme

Before introducing our improved memory cutoff, we first review the ADT scheme and the approximations required. This will also aid us in introducing our notation. The class of models to which the scheme is applicable are those consisting of a small system of interest linearly coupled to a macroscopic bath of bosonic modes. The generic Hamiltonian of such models is

[TABLE]

where $H_{0}$ is the free system Hamiltonian and $H_{B}$ contains both the free bath Hamiltonian and the system-bath interaction. Here, $a_{\alpha}^{\dagger}$ ($a_{\alpha}$) and $\omega_{\alpha}$ are the creation (annihilation) operators and frequencies of the $\alpha$th oscillator. The system operator $\hat{s}$ couples to the bath operators $\hat{B}_{\alpha}=g_{\alpha}a_{\alpha}+g_{\alpha}^{*}a_{\alpha}^{\dagger}$ with coupling constants $g_{\alpha}$.

To simplify our notation we work in the Liouvillian representation such that operators in Hilbert space are represented by vectors in Liouville space. To parameterise the $D$ dimensional Hilbert space of the system we use the $D$ eigenstates of $\hat{s}$. Operators in this space are vectorized in the following way

[TABLE]

where the sum over $S$ runs over the $D^{2}$ pairs of $\{s^{+},s^{-}\}$ and $\hat{s}\ket{s^{+}}=s^{+}\ket{s^{+}}$ and likewise for $s^{-}$. We use the notation $\left.\Ket{x}\right\rangle$ to mean a vector in Liouville space. The bath Hilbert space can also be represented in a similar way, though we do not need to define its basis explicitly in what follows. The evolution of the reduced system, assuming factorizing initial conditions is now represented as:

[TABLE]

with the Liouvillian $\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{B}$, where $\mathcal{L}_{0}$ and $\mathcal{L}_{B}$ generate coherent evolution caused by $H_{0}$ and $H_{B}$ respectively. Recently it has been shown Barth et al. (2016) that additional Markovian non-unitary dynamics of the reduced system can be incorporated by adding Lindblad-type dissipators into the free system Louivillian $\mathcal{L}_{0}$, making it straightforward to account for coupling to other baths for which the Born-Markov approximations are well justified. In addition to factorising initial conditions we also assume the initial state of the bath is that of thermal equilibrium when no system is present $\rho_{B}=\exp(-\sum_{\alpha}\omega_{\alpha}a^{\dagger}_{\alpha}a_{\alpha}/T)/\mathcal{Z}$, with temperature $T$ and partition function $\mathcal{Z}$.

The first approximation made to make Eq. (4) computable is to factorize the long time propagator into $N$ short time propagators $\text{e}^{\mathcal{L}t}=(\text{e}^{\mathcal{L}\Delta t})^{N}$ and then to employ a Trotter splitting between the system and bath parts Trotter (1959)

[TABLE]

on each of these. The error introduced in this process is $\mathcal{O}(\Delta t^{2})$. We note that the argument that now follows can be easily adapted to use a symmetrized Trotter splitting Suzuki (1976); Makri and Makarov (1995a, b) that improves the error to $\mathcal{O}(\Delta t^{3})$. All the numerical results we present do include this symmetrized splitting, but for simplicity of notation we will use the definition in Eq. 5 here.

Tracing out the bath degrees of freedom then results in a reduced density matrix at time $t_{N}=N\Delta t$, whose elements are

[TABLE]

The functions $F\left(\{S_{k}\}\right)$ and $I\left(\{S_{k}\}\right)$ constitute the free part of the evolution and the discretized Feynman influence functional respectively, and are given by

[TABLE]

The coefficients $\eta_{k-k^{\prime}}$ quantify the non-Markovian ‘interaction’ between the reduced system at different times $t_{k}$ and $t_{k^{\prime}}$ and are defined as

[TABLE]

in terms of the bath autocorrelation function

[TABLE]

where $J(\omega)=\sum_{\alpha}|g_{\alpha}|^{2}\delta(\omega_{\alpha}-\omega)$ is the spectral density of the bath.

In its current form Eq. (6) is, in principle, numerically computable, though would require the storage of the $D^{2N}$ numbers making up the discretized influence functional. The exponential dependence of storage requirements on the number of timesteps places strict limits on the length of time for which a given simulation can be propagated. To circumvent this issue we make use of the finite memory approximation originally introduced by Makri and Makarov (1995a, b). For a continuum of bath oscillator modes the correlation function Eq. (11) decays at worst algebraically and in finite time becomes effectively zero. This in turn implies that the coefficients $\eta_{k-k^{\prime}}$ decay to zero as the difference $t_{k}-t_{k^{\prime}}$ is increased. The finite memory approximation involves setting those coefficients for which $k-k^{\prime}>\Delta k$ such that the time difference $t_{k}-t_{k^{\prime}}>\Delta k\Delta t\equiv\tau_{c}$ identically equal to zero. Here we have introduced $\tau_{c}$, the bath cutoff time, which must be suitably large to ensure converged results. In Fig. 1 we visualize this memory cutoff in terms of two dimensional integral regions over $C(t^{\prime}-t^{\prime\prime})$. For fixed $k$ the coefficients $\eta_{k-k^{\prime}}$ correspond to columns of a number $\Delta k$ square integral regions extending up to the triangular regions which give $\eta_{0}$. The finite memory approximation is then made by setting all of the coefficients in the green region to zero. This has the effect that the dynamics we calculate in region A are exact and the approximation only effects what happens for times $t>\tau_{c}+\Delta t$.

With the finite memory approximation we can now reformulate Eq. (6) in terms of the iterative propagation of an object known as the augmented density tensor (ADT as before). This propagation is computationally realised most efficiently as contractions between multidimensional tensors but the analogy to standard density matrix propagation is most easily seen by representing it as operator-vector multiplication in an augmented Liouville space:

[TABLE]

where $\left.\Ket{\rho_{R}^{A}}\right\rangle$ is the vectorized ADT and $\mathcal{L}^{A}$ is the augmented Liouvillian. The augmented space is constructed by taking the product of $\Delta k$ copies of the original Hilbert space. The basis of this space at a given time $t_{k}$ is constructed by taking a product of reduced system basis states at that time and the previous $\Delta k-1$ times $t_{k-1}\ldots t_{k-\Delta k+1}$

[TABLE]

thus giving the augmented Hilbert space a dimension of $D^{2\Delta k}$. The information stored in the ADT is a set of amplitudes weighting each of the trajectories the reduced system could have taken through its Hilbert space in the previous $\Delta k$ timesteps of the evolution. The physical reduced density matrix at a given time is then found by summing these amplitudes:

[TABLE]

Computing Eq. (6) with the memory cutoff requires the storage of the $D^{2\Delta k}$ elements of the vector representing the augmented state and the $D^{4\Delta k}$ elements of the propagator matrix. As mentioned above this is not the optimal representation for actually carrying out the propagation; in fact the alternative tensor representation of the propagator only has $D^{2(\Delta k+1)}$ elements. In Appendix A we explicitly construct this more efficient tensor propagation and also give the matrix elements of the propagator, $\left\langle\Bra{S^{A}_{k}}\right.e^{\mathcal{L}_{A}\tau_{c}}\left.\Ket{S^{A}_{k-\Delta k}}\right\rangle$, as well as the components of the initial augmented state vector, $\left\langle\Braket{S^{A}_{\Delta k}}{\rho_{R}^{A}(\tau_{c})}\right\rangle$, required for both tensor and matrix representations of the method.

The scheme then has two parameters which need to be adjusted to ensure convergence of results. The memory cutoff time, $\tau_{c}$, should be made large enough, by increasing $\Delta k$, such that increasing it any further has negligible effects on the final result. At the same time the error induced by the Trotter splitting must also be eliminated by decreasing the timestep $\Delta t$, which in turn decreases $\tau_{c}$ for a given $\Delta k$. Thus, achieving the best results requires both maintaining a small enough timestep $\Delta t$ to eliminate the Trotter error while at the same time keeping it large enough such that the number of timesteps required to capture the correlation time of the bath is kept as small as possible. Details on the method for finding the optimal value of $\Delta t$ that minimises the overall error can be found in Ref. Eckel et al. (2006).

Finally, we note that the ADT has many properties in common with a standard density matrix: it is both Hermitian and has unit trace, from which it follows that a system density matrix calculated from it is also guaranteed to have these properties. To describe a physical state it is also necessary that the system density matrix is positive. Proving this is, in general, more difficult and it is known that the ADT scheme does not guarantee positivity of the reduced system state Vagov et al. (2011). In the following section we will gain an understanding of how non-positive reduced system states can occur and describe a procedure to help prevent this from happening.

III Memory Cutoff in an Exactly Soluble Model

To investigate the role of the finite memory approximation in the occurrence of non-positive reduced system states, here we study the effects of making this approximation on an exactly soluble pure dephasing model where $[H_{0},\hat{s}]=0$. The simplest example of such a model is that of a two level system described by the Hamiltonian $H_{0}=\epsilon\sigma_{z}/2$ which couples to the bath via $\hat{s}=\sigma_{z}$. This is the independent boson model Breuer and Petruccione (2002) which we will use as an example here. Note however that what follows can be applied to any model which satisfies the commutation relation above.

III.1 Independent Boson Model

For this model the Trotter splitting in Eq. (5) is exact since the system and bath Hamiltonians commute. Also, by moving to the interaction picture of the reduced system we may ignore the dynamics induced by $H_{0}$ and the free system propagator is therefore the identity

[TABLE]

The summation in Eq. (6) can now be carried out exactly (without the finite memory approximation) to obtain

[TABLE]

where

[TABLE]

Here we have defined the function

[TABLE]

which governs all dynamics induced by the interaction with the bath. We see that state populations on the diagonal of $\rho_{R}(t)$ (where $s^{+}=s^{-}$) do not undergo evolution and all dynamics are in the decay and oscillations of the coherences.

In what follows we examine the effect of making the finite memory approximation on the function $\eta(t)$. Note that this also describes the effects of making this approximation on the discretized influence functional since the $\eta_{k-k^{\prime}}$ coefficients can be written entirely in terms of this function,

[TABLE]

To work out how the solution is changed when we impose the memory cutoff we must consider what happens to the function $\eta(t)$ when we remove those coefficients for which $k-k^{\prime}>\Delta k$ from the sum in Eq. (19). In Fig. 1 we see there are two distinct regions of the integral domain left after performing the memory cutoff: when $t<\tau_{c}+\Delta t$, (region A) $\eta$ is unaffected by the cutoff and region B where $t>\tau_{c}+\Delta t$ and at least some correlations are cut off. The sum of the coefficients within region A is exactly $\eta(\tau_{c}+\Delta t)$, while using Eq. (20) the sum of all the coefficients in a single strip of region B is $\eta(\tau_{c}+\Delta t)-\eta(\tau_{c})$. There are $N-\Delta k-1$ of these strips in total and so by summing up all coefficients remaining after the finite memory approximation we find the function $\eta(t)$ is approximated as

[TABLE]

This result can also be obtained in the $\Delta t\to 0$ limit by setting $C(t)=0$ for $t>\tau_{c}$ in the integral of Eq. (19). By making a change of variable $s=t^{\prime}-t^{\prime\prime}$, we find

[TABLE]

where the dot indicates a time derivative. Thus the dynamics of the pure dephasing model with the finite memory approximation imposed are as follows: At times $t\leq\tau_{c}+\Delta t$ the exact dynamics are followed. While for $t>\tau_{c}+\Delta t$ there is an exponential decay with rate $\gamma_{c}=\operatorname{\mathrm{Re}}[\dot{\eta}(\tau_{c})]$ accompanied by a Lamb-shift of the system frequencies, $\Lambda_{c}=\operatorname{\mathrm{Im}}[\dot{\eta}(\tau_{c})]$.

It is easy to see now how a problem can arise if the memory cutoff is at a point where the gradient of $\eta(t)$ is negative. The resultant decay rate is negative and hence there is exponential growth of all coherences for $t>\tau_{c}$. This inevitably leads to the system density matrix becoming non-positive.

Additionally, unless $\operatorname{\mathrm{Re}}[\dot{\eta}(\tau_{c})]=0$ identically, finite steady state coherences will be impossible to produce. Increasing $\tau_{c}$ increases the time it takes for the coherences to decay or grow away from the true steady state, but for finite $\tau_{c}$ the convergence of the finite memory approximation, and hence the ADT scheme, cannot be guaranteed at arbitrary times.

To investigate how these problems manifest in the dynamics of a real example we need to specify the environmental coupling, quantified by the spectral density $J(\omega)$, which determines the form of $\eta(t)$. We consider the spectral density

[TABLE]

which is characterised by the dimensionless coupling constant, $\alpha$, ohmicity, $\nu$ and an exponential cutoff at a scale given by $\omega_{c}$. It has been found that the general condition for $\dot{\eta}(t)>0$ is that the function $J(\omega)\coth(\omega/2T)/\omega^{2}$ be convex, and that for spectral densities of the form Eq. 25 this condition is fulfilled when $\nu<\nu_{\text{crit}}$ for some critical ohmicity $\nu_{\text{crit}}$ Haikka et al. (2013). The value of $\nu_{\text{crit}}$ varies monotonically from $\nu_{\text{crit}}=2$ and $\nu_{\text{crit}}=3$ as temperature is increased from 0 to $\infty$. Negative decay rates arise when $\nu>\nu_{\text{crit}}$. Therefore at $T=0$ it is in the superohmic regime, where $\nu>2$, that we expect to find the finite memory cutoff results in non-positive dynamics. To demonstrate this we plot in Fig. 2 the time dependent decay rate resulting from both ohmic and superohmic spectral densities at zero temperature. We also plot the resultant dynamics of the coherence, for both the exact solution of the independent boson model, Eq. (17), and with the memory cutoff imposed. In this case the imaginary part of the exponential in Eq. (17) disappears, since the eigenvalues of $\hat{s}$ have equal magnitude, so the dynamics are that of pure decay. For the ohmic case the decay rate is always positive and using the memory cutoff produces exponential decay to the correct steady-state of zero coherence. There is still a significant discrepancy at intermediate times, but this can be reduced by increasing $\tau_{c}$. For a superohmic spectral density there is a large window of time when the decay is negative. Making the cutoff in this window results in unphysical dynamics with exponential growth. Notably, in this superohmic case the decay rate approaches zero from below asymptotically. This means that increasing $\tau_{c}$ does not remove this spurious asymptotic behaviour it just shifts its onset to longer and longer times.

To conclude this analysis we point out a situation in which even the ohmic and sub-ohmic regimes $0<\nu\leq 2$ can display negative decay rates and hence non-positive states. In constructing physically realistic models one may want to account for how spatially separated states interact with the same bath of oscillators. For the case of the two-level independent boson model considered above in three dimensions the spatial separation of the two system eigenstates is accounted for via a multiplicative correction to the spectral density $\propto 1-\operatorname{\mathrm{sinc}}(k\omega)$ McCutcheon and Nazir (2011). Here $k=d/c$ where $d$ is the distance between the two sites and $c$ is the speed at which the bath excitations propagate. Thus, for small separation distance, $k\omega_{c}<1$, spectral densities of the form Eq. (25) go as $J(\omega)\propto\omega^{\nu+2}e^{-\omega/\omega_{c}}$ and the effective critical ohmicity for the transition to non-Markovian dynamics with negative decay rates now lies in the range $0<\nu_{crit}\leq 1$ and the potential for the memory cutoff to produce unphysical dynamics exists for all values of $\nu$.

III.2 Fixing the non-positive evolution

We now seek a way of taking the memory cutoff such that the exact solution Eq. 17 is always reproduced using the ADT method. A general way to do this would be to define a new set of the coefficients $\tilde{\eta}_{k-k^{\prime}}$, $k-k^{\prime}\leq\Delta k$ in such a way that their sum is constrained as follows,

[TABLE]

thus reproducing the exact dynamics governed by $\eta(t)$, independent of both $\tau_{c}$ and $\Delta t$. It seems reasonable that we should attempt to redefine as few of the coefficients as possible, since they already exactly account for non-Markovian correlations within a timespan $\tau_{c}$ (up to the Trotter error). The key idea behind the finite memory cutoff is that the most non-local temporal correlations in the system are the most insignificant, hence $|\eta_{k-k^{\prime}}|\approx 0$ for large $k-k^{\prime}$. Thus we redefine only the coefficients with the largest $k-k^{\prime}=\Delta k$ as follows

[TABLE]

This redefinition will not introduce any problems as long as $|\eta_{k-k^{\prime}}|\approx|\tilde{\eta}_{k-k^{\prime}}|$ which is true since $\tau_{c}$ is already large enough to naïvely justify the finite memory approximation. This redefinition is visualized in Fig. 3 in a similar way to the standard memory cutoff in Fig. 1. From this schematic one can identify that this redefinition corresponds to extending the lower limit of the integral down to zero. This means that it is straightforward to calculate the redefined coefficient:

[TABLE]

In terms of implementing the ADT scheme the main consequence is that now a new propagator must be constructed at every step of the iterative propagation, though the actual structure of each propagator is essentially identical to the original. Thus, although slightly more time consuming, carrying out the propagation is no more complicated than before. The time dependence in our method is fundamentally different from that generated by time-dependent Hamiltonians, in that it is due to the non-local influence of the bath interaction itself. In the standard finite memory approximation the finite ranged non-Markovian correlations allow the problem to be reformulated as one with Markovian evolution but in a higher dimensional augmented Hilbert space as in Eq. (13). We have gone one step further by allowing $\mathcal{L}_{A}$ to be time dependent, thus allowing transient non-positive evolution of the reduced system state beyond the cutoff time $\tau_{c}$ while still maintaining overall positive evolution from the initial state. This is not possible in the standard ADT scheme where the augmented Liouvillian is time independent and so must have no positive eigenvalues to ensure positive evolution. We point out here that the instability of steady states and non-positivity in the ADT scheme for superohmic environments has been recognized before, and a similar redefinition of the $\eta_{k-k^{\prime}}$ coefficients has been proposed Vagov et al. (2011); G Nahri1 et al. (2016) but only in a way that gives qualitatively correct asymptotic behaviour of the pure dephasing model above, but does not reproduce the exact solution identically.

In the following section we show that this new method is useful for more complicated models which do not have an exact solution by applying it to the more general spin-boson model.

IV Application to the Spin-Boson Model

The spin-boson model Leggett et al. (1987) provides a paradigmatic example of an open quantum system and can be used to model a wide variety of physical systems. For example, it has been applied to the problem of electron transfer in biological aggregates Xu and Schulten (1994), exciton dynamics in quantum dots Ramsay et al. (2010a), transport in mesoscopic systems Brandes (2005) and chemical reactions Garg et al. (1985). The system Hamiltonian is $H_{0}=\epsilon\sigma_{z}+V\sigma_{x}$ and the coupling to the environment is through the operator $\hat{s}=\sigma_{z}$. Here $\epsilon$ again gives the energy splitting of the two levels, while the $\sigma_{x}$ term generates coherent transitions between the two states. The addition of this $\sigma_{x}$ term breaks the integrability of the independent boson model and allows for a rich variety of physics to be explored.

In what follows we will be particularly interested in the case where the environment is superohmic, since this is where we found the most pathological behaviour in the independent boson model. The superohmic regime is most studied in the context of quantum dots strongly coupled to a phononic environment Ramsay et al. (2010b, a); Kok and Lovett (2010). For many parameters a polaron master equation provides a successful route to capturing non-perturbative effects Nazir and Mccutcheon (2016), but for highly non-Markovian environments this approach fails. Here we show how the ADT scheme is able to capture backflow of energy from the environment to the system in this regime.

To show the improvement gained by the new integration scheme which we detailed in the previous section, we show how the convergence of the results with $\Delta k$ changes as compared to the standard approach. In Fig. 4 we study the spin-boson model with $\epsilon=0$ at zero temperature with a large reservoir cutoff frequency $\omega_{c}/V=5$. To find an approximation for the error in our results we compare everything to the most converged case: that using the new memory cutoff at $\Delta k=14$. We then calculate the time average of the trace distance Breuer et al. (2016) between these converged results and each other case. The trace distance is given by

[TABLE]

where where $|A|=\sqrt{AA^{\dagger}}$, $\rho_{1}(t)$ is the converged density matrix and $\rho_{2}(t)$ is the density matrix we want to compare to $\rho_{1}(t)$. We show this deviation as a function of $\Delta k$ for both approaches in the panel (c) of Fig. 4. We see that our new approach converges much more quickly than the standard finite memory approximation: by $\Delta k=11$ the errors in the new approach are negligible while the standard finite memory approximation still has significant errors at $\Delta k=14$. This can also be seen in the dynamics in Fig 4(a) and (b) where we see that our new approach converges much more quickly when increasing $\Delta k$.

It is difficult to find convergence with this set of parameters because the gradient of $\eta(t)$ approaches zero from below (as can be seen in the inset to Fig. 4(c)) and so the problems we found for the independent boson model still occur, although here they are less severe and the standard approach only gives unphysical results at very small $\Delta k$. This problem becomes even worse if we move further towards the so-called scaling limit of the model Leggett et al. (1987), where $\omega_{c}$ is the largest energy scale in the problem, by increasing the value of $\omega_{c}$. This means a smaller required $\tau_{c}$, but also results in a larger magnitude of $\dot{\eta}(\tau_{c})$ so that the timescale for the onset of divergent dynamics becomes much shorter and achieving convergence over appreciable lengths of time becomes very difficult.

Now we have shown the improvement gained by our improved memory cutoff we examine in detail the dynamics of the spin-boson model in a parameter regime which is difficult to access using the standard memory cutoff approximation: when there is non-Markovian behaviour far from the scaling limit with a lower cutoff frequency $\omega_{c}/V=1$ and at finite temperature. In this limit the required value of $\Delta k$ without our improvement would be much too large to store in memory for a standard computer.

In Fig. 5(a) we show the population difference dynamics with initially excited spin, $\rho_{0}=\ket{1}\bra{1}$, for system-bath couplings between $0.5<\alpha<1.9$. At low coupling strengths we see simple damped oscillations, as would be captured by a weak coupling master equation. Increasing the coupling changes these oscillations from underdamped to overdamped, as would be found using for example a polaron master equation. However, the overdamped regime is accompanied by a highly non-Markovian feature: a revival of the population difference which becomes more pronounced at stronger couplings.

In order to establish a more concrete quantification of the non-Markovian nature of these dynamics we examine how distinguishable two distinct initial states are as they evolve in time. To do this we again use the trace distance $D(\rho_{1}(t),\rho_{2}(t))$ as a natural measure. For a Markovian system as time evolves all states move closer to the steady state and so asymptotically we expect $\lim_{t\to\infty}\rho_{1}(t)=\lim_{t\to\infty}\rho_{2}(t)=\rho_{\text{steady}}$, provided the steady state $\rho_{\text{steady}}$ is unique. This means that $D(\rho_{1}(t),\rho_{2}(t))$ monotonically decays to 0 for any two initial states. If, however, the dynamics are non-Markovian then, for certain pairs of initial states, there can be increases in the trace distance as a function of time as the environment allows information to flow back into the system. It is therefore evident that quantity of interest is actually

[TABLE]

which is positive when these non-Markovian features occur.

In Fig. 5(b) and (c) we show the trace distance and $\Delta I(t)$, for the same set of parameters as in (a), using the orthogonal initial states $\rho_{1}(t)=\ket{1}\bra{1}$ and $\rho_{2}(t)=\ket{-1}\bra{-1}$. At low values of coupling there are weak oscillations in $\Delta I(t)$ between negative and positive which are damped and pushed out to longer times as the coupling is increased. Increasing the coupling still further leads to an additional positive region of $\Delta I(t)$ at short times, which remains there, closely following the early time revival we see in the population differences in Fig. 5. This observation confirms that this feature is a signature of information backflow. Similar features have been observed in biased spin-boson models (i.e. $\epsilon\neq 0$) Prior et al. (2010); Wilner et al. (2015); Jang (2009). Here, we show quantitatively that non-Markovian features are present in the dynamics of the unbiased spin-boson model.

To analyse the dependence of the non-Markovian revival in the population dynamics on the parameters describing the bath we fit the dynamics to the sum of two decaying oscillations:

[TABLE]

with the constraint $\alpha_{z}(0)=A+B\sin(\phi)=1$. The frequency $\omega_{1}$ and decay rate $\gamma_{1}$ parameters capture the weak-coupling and long time dynamics, while $\omega_{2}$ and $\gamma_{2}$ describe the higher frequency short time dynamics of the revival. An example of one of these fits, for $\alpha=0.6$, is shown in the inset to Fig. 6.

We show how $\gamma_{2}/V$ and $\omega_{2}/V$ vary as the coupling to the bath, $\alpha$, is increased in Fig. 6(a). Both $\gamma_{2}/V$ and $\omega_{2}/V$ increase with coupling at roughly the same rate and, although the behaviour is certainly not linear, this dependence implies that increasing coupling simply causes the timescale over which the revival occurs to shorten. This is consistent with the revival being due to a backflow of information from the bath: we would expect larger couplings to simultaneously cause quicker information backflow and to wash it out more quickly. In Fig. 6(b) we show what happens as the temperature of the bath is changed. Increasing temperature both increases the effective coupling to the bath but also scrambles temporal correlations, reducing the bath correlation time. Therefore we see both the damping and frequency of the revival increase with increasing temperature.

Before concluding, we compare the results above to those obtained using a Markovian polaron master equation Nazir and Mccutcheon (2016). If our interpretation of the revival is correct, then it is obvious that any approach which does not fully account for the dynamics of the bath will be unable to reproduce this feature. For example, the population difference dynamics predicted using the polaron master equation are of the same form as Eq. (31) but with $\gamma_{1}=\gamma_{2}$ and $\omega_{1}=\omega_{2}$ McCutcheon and Nazir (2011); Nazir and Mccutcheon (2016). This is shown in Fig. 7. For the largest values of $\omega_{c}$ in (a) and (b) the system is in the scaling limit Leggett et al. (1987), where the polaron technique works well, and indeed we find good agreement of the two methods, with only small deviations occurring at very short timescales. As $\omega_{c}$ decreases, we enter a more non-Markovian regime and we see that the polaron master equation starts to fail. It is difficult to find converged results which show this effect more clearly than shown here, since in the scaling limit divergent dynamics occur on very short timescales, as discussed earlier. The long time dynamics predicted by the polaron equation are qualitatively correct, taking the form of an exponentially decaying oscillation, but as anticipated it completely fails to capture the non-Markovian revival which we found above using the improved ADT. We attribute this to the breakdown of the Markov approximation used in deriving the polaron master equation. From our earlier discussion, we know that the revival in dynamics is most prominent for smaller values of $\omega_{c}$, i.e. for baths with longer correlation times. This, together with the poor agreement of ADT with the polaron master equation at short times, confirms our conclusion that this feature is a highly non-Markovian effect which could not be produced using any form of time-local master equation.

V Summary

We have shown how the standard finite memory approximation in the ADT numerical scheme can cause unphysical behaviour resulting in periods of non-positive evolution. We have provided an improvement to this method which is able to reproduce the exact solution to the pure dephasng model. To demonstrate the applicability of this new method we have shown converged results for dynamics of the symmetric spin-boson model in a superohmic environment. We found that our new method is able to reach convergence using significantly smaller computational resources than the standard finite memory approximation. At strong system-bath coupling strengths in this regime we found highly non-Markovian revivals in the dynamics which are accompanied by non-monotonicity of the trace distance between different initial states. These features are not present in simpler analytical approaches such as polaron master equations.

Acknowledgements

We acknowledge useful discussions with E. M. Gauger. We thank D. P. S. McCutcheon for providing code used for the polaron simulations. AS acknowledges a studentship from EPSRC (EP/L505079/1). BWL acknowledges support from EPSRC (EP/K025562/1). PGK acknowledges support from EPSRC (EP/M010910/1). The research data supporting this publication can be found at http://dx.doi.org/10.17630/21764101-6493-46e7-ab85-65cccc1fb9e5.

Appendix A Explicit Construction of Propagation Methods

We start by rewriting the summand of Eq. (6) in the following form

[TABLE]

where

[TABLE]

with

[TABLE]

There is no $S_{0}$ dependence in any of the $\tilde{I}(S_{k},S_{k^{\prime}})$ so the $S_{0}$ summation in Eq. (6) can be carried out, propagating the reduced system for a time $\Delta t$ under the free Liouvillian $\mathcal{L}_{0}$ into the state $\rho_{R}(\Delta t)$. The finite memory approximation, setting $\eta_{k-k^{\prime}}=0$ for $k-k^{\prime}>\Delta k$, then translates to $\tilde{I}(S_{k},S_{k^{\prime}})=1$ for $k-k^{\prime}>\Delta k$. This allows the product over $\tilde{I}(S_{k},S_{k^{\prime}})$’s to be rewritten as

[TABLE]

where we have defined

[TABLE]

to be the elements of the propagator tensor, and

[TABLE]

to be the elements of the initial augmented density tensor. The summations over the $S_{k}$ can now be carried out one at a time by iteratively multiplying and contracting tensors:

[TABLE]

with the initial condition Eq. (37). The reduced system density matrix at a time $t_{k}$ is then retrieved by summing over all but the $S_{k}$ index

[TABLE]

The propagator tensor Eq. (36) has $\Delta k+1$ indices and so has $D^{2(\Delta k+1)}$ elements, making the tensor contraction representation of the propagation less demanding than representing the augmented state as a vector and propagating with a $D^{4\Delta k}$ sized matrix. Within the Liouvillian matrix formalism for propagation, Eq. (13), and using the basis defined in Eq. (14), the elements of the augmented state at a time $t_{k}$ are:

[TABLE]

and the matrix elements of the propagator across the timespan $\tau_{c}$ have the simple analytic form:

[TABLE]

Thus another way of solving the problem would be to construct and diagonalize this propagator to find the eigenvectors of the augmented Liouvillian $\mathcal{L}^{A}$, though for typical values $\Delta k\sim 10$ it is much more efficient to use the iterative propagation schemes above.