A hybrid memory kernel approach for condensed phase non-adiabatic dynamics

TL;DR

This paper introduces a hybrid memory kernel method for non-adiabatic dynamics in condensed phase systems, combining exact slow-bath solutions with NIBA to improve accuracy in challenging regimes.

Contribution

It develops a hybrid approach that enhances second-order methods with higher-order effects, improving dynamics predictions for complex systems.

Findings

The hybrid kernel outperforms NIBA in large diabatic coupling regimes.

The method maintains computational efficiency similar to NIBA.

It provides an exact solution in the slow-bath limit and effective approximations for faster baths.

Abstract

The spin-boson model is a simplified Hamiltonian often used to study non-adiabatic dynamics in large condensed phase systems, even though it has not been solved in a fully analytic fashion. Herein, we present an exact analytic expression for the dynamics of the spin-boson model in the infinitely slow bath limit and generalize it to approximate dynamics for faster baths. We achieve the latter by developing a hybrid approach that combines the exact slow-bath result with the popular NIBA method to generate a memory kernel that is formally exact to second order in the diabatic coupling but also contains higher-order contributions approximated from the second order term alone. This kernel has the same computational complexity as NIBA, but is found to yield dramatically superior dynamics in regimes where NIBA breaks down---such as systems with large diabatic coupling or energy bias. This indicates that this hybrid approach could be used to cheaply incorporate higher order effects into second order methods, and could potentially be generalized to develop alternate kernel resummation schemes.

Equations

$${\mathbf{H}}_{SB}$$

$${\mathbf{h}}_1$$

$${\mathbf{h}}_2$$

$$J(\omega )=\sum _{i=1}^N\frac{\pi c_i^2}{2m_i{\omega }_i}\delta (\omega -{\omega }_i)$$

$$\frac{d{p}_1}{dt}$$

$$\frac{d{p}_2}{dt}$$

$$\frac{d⟨{\sigma }_z⟩}{dt}$$

$$⟨{\sigma }_z⟩(t)$$

$$K_{\pm }(t,t_1)$$

$$f_{2n}(t_1,t_2,t_3\dots t_{2n})$$

$$O(t)$$

$$\frac{d}{dt}{⟨{\sigma }_z⟩}^{(2)}(t)$$

$$K_{11}^{(2)}(t,t_1)$$

$$K_{22}^{(2)}(t,t_1)$$

$$Q^{\prime }(t)$$

$$Q^{\prime \prime }(t)$$

$$K_{11}^{(2)}(t,t_1)$$

$$K_{22}^{(2)}(t,t_1)$$

$$\varphi (t,t_1)$$

$$\ln f_{2n}(t_1,t_2,t_3\dots t_{2n})$$

$$f_4(t_1,\dots t_4)$$

$$\ln f_2(t_1,t_2)$$

$$b$$

$$a$$

$$f_2(t_1,t_2)$$

$$f_{2n}(t_1,t_2\dots t_{2n})$$

$$\dot{⟨{\sigma }_z⟩}(t)=V^2{\dot{⟨{\sigma }_z⟩}}^{(2)}(t)-2V^3\underset{0}{\overset{t}{\int }}{J}_1(2V{t}_1){\dot{⟨{\sigma }_z⟩}}^{(2)}\left(\sqrt{t^2-t_1^2}\right)d{t}_1$$

$${\dot{⟨{\sigma }_z⟩}}^{(2)}(t)=-\mathbf{Re}\left[\frac{2\sqrt{\pi }{e}^{-\frac{b^2}{4a}}\text{erf}\left(\frac{2at+ib}{2\sqrt{a}}\right)}{\sqrt{a}}\right]$$

$$K_{11/22}^I(t,t_1)$$

$$+\left(1-u(t-t_1)\right)K_{11/22}^{(2)}(t,t_1)$$

$$\dot{⟨{\sigma }_z⟩}(t)=$$

$$-2V^3\underset{0}{\overset{t}{\int }}{J}_1(2Vf(t_1)t_1){\dot{⟨{\sigma }_z⟩}}^{(2)}\left(\sqrt{t^2-t_1^2}\right)d{t}_1$$

$${\mathbf{\rho }}_I(t)$$

$$={\mathbf{\rho }}_I(0)-iV\underset{0}{\overset{t}{\int }}[{\sigma }_I^x(t_1),{\mathbf{\rho }}_I(0)]d{t}_1-V^2\underset{0}{\overset{t}{\int }}d{t}_1\underset{0}{\overset{t_1}{\int }}[{\sigma }_I^x(t_1),[{\sigma }_I^x(t_2),{\mathbf{\rho }}_I(0)]]d{t}_2\dots$$

$$=\sum _{m=0}^{\infty }{V}^m{\mathbf{\rho }}_I^{(m)}(t)$$

$$p_1(t)$$

$$=\sum _{m=0}^{\infty }{V}^m{\mathbf{Tr}}_B\left[{\mathbf{\rho }}_I^{(m)}(t)|1⟩⟨1|\right]$$

$$p_2(t)$$

$$f_2(t_1,t_2)$$

$$\frac{\partial f_2(t_1,t_2)}{\partial t_1}$$

Full text

I Introduction

Condensed phase chemical processes are strongly influenced by the large number of solvent degrees of freedom that interact with the reacting moleculesMarcus (1956, 1965, 1993); Barbara, Meyer, and Ratner (1996). Accurate modeling of condensed phase dynamics thus requires incorporation of solvent effects, but this is difficult to achieve via atomistic ab initio techniques as the computational cost scales steeply with the number of solvent coordinates treated quantum mechanically. The generality of such atomistic approachesTully (1990); Cao and Voth (1994a, b); Ben-Nun and Martınez (1998); Ben-Nun, Quenneville, and Martínez (2000); Craig and Manolopoulos (2004) makes them extremely useful nonethelessHammes-Schiffer and Tully (1994); Toniolo, Granucci, and Martínez (2003); Toniolo et al. (2004); Menzeleev, Ananth, and Miller III (2011), but they may not represent the best tools for many problems.

The success of Marcus theoryMarcus (1956, 1965, 1993), on the other hand, suggests that simple model systems containing a few parameters (obtained from experiment or calculations on single potential energy surfaces) could be an alternate route to study condensed phase non-adiabatic dynamics. A well known example of such a model is the spin-boson modelCaldeira and Leggett (1983), which is used quite often to study the dynamics of electron transferBarbara, Meyer, and Ratner (1996); Bader, Kuharski, and Chandler (1990); Xu and Schulten (1994); Evans, Nitzan, and Ratner (1998) as it reduces to the highly successful Marcus theory under appropriate limitsXu and Schulten (1994). This suggests that it can be viewed as a dynamical generalization of Marcus theory that could be useful for studying dynamics of electron transfer and non-adiabatic molecular dynamics in general.

However, an exact analytic solution for the population dynamics of the spin-boson model has not yet been found, despite the relative simplicity of the Hamiltonian. Various numerical approaches like Quasi-Adiabatic Propagator Path Integral (QUAPI) Makarov and Makri (1993, 1994); Topaler and Makri (1994); Makri (1995); Makri and Makarov (1995a, b), Multi-Configurational Time-Dependent Hartree (MC-TDH)Wang, Thoss, and Miller (2001); Thoss, Wang, and Miller (2001), Hierarchical Equations of Motion (HEOM)Ishizaki and Tanimura (2005); Tanimura (2006) and Quantum Monte Carlo (QMC)Egger and Weiss (1992); Winter et al. (2009) have therefore been used to study the dynamics of the spin-boson model. Many of these approaches are formally numerically exact, but often incur significant computational cost that makes exploration of parameter space difficult. Analytic approximations like the Non-Interacting Blip Approximation (NIBA)Caldeira and Leggett (1983); Chakravarty and Leggett (1984); Aslangul, Pottier, and Saint-James (1985) and generalized quantum master equationsSparpaglione and Mukamel (1988); Mavros and Van Voorhis (2014); Kelly et al. (2016) are also often used to study the dynamics of the spin-boson model, although they may also be quite computationally expensive or be effective only in certain parameter regimes.

In this paper, we present an exact analytic solution to the spin-boson model in the limit of an infinitely slow bath, and combine this result with NIBA to generate a hybrid kernel for the dynamics of general baths with arbitrary spectral densities. This hybrid kernel approximation has the same computational complexity as NIBA, but is found to yield much superior results in regimes where the latter fails, being approximately as effective as a fourth order resummed memory kernel approach for all cases tested. In particular, it is found to give quite good performance in problematic cases involving low temperatures, large diabatic coupling and large energy biases, where NIBA is known to give qualitatively incorrect behavior. This indicates that the hybrid approach is a promising alternative to naive NIBA, as it has the same computational complexity but yields superior results in all regimes tested. Consequently, it can be used to cheaply explore a large parameter space for screening purposes over a more expensive method like MC-TDH. The developed approach is also able to approximate dynamics from second order information alone, suggesting that it can potentially be generalized to develop a suite of resummation methods that can be applied to problems more general than the spin-boson model.

II The spin-boson model

The spin-boson modelCaldeira and Leggett (1983) is one of the simplest models for open quantum systems, describing the interaction between a two level system and a bath of harmonic oscillators. The model assumes that the diabatic coupling $V$ between the two system states is constant and the system-bath coupling is linear in the bath coordinates. There may however be a non-zero intrinsic energy bias $\epsilon$ between the two system states.

Therefore, the Hamiltonian $\mathbf{H}_{SB}$, in the basis of orthogonal diabatic states $\left|1\right\rangle=\begin{pmatrix}1\\ 0\end{pmatrix}$ and $\left|2\right\rangle=\begin{pmatrix}0\\ 1\end{pmatrix}$, is given by:

[TABLE]

The harmonic bath has a set of modes $\{x_{i}\}$ with frequencies $\{\omega_{i}\}$ and masses $\{m_{i}\}$. These modes interact linearly with the system states $\left|1\right\rangle$ and $\left|2\right\rangle$ via the coupling constants $\{c_{i}\}$. Consequently, $\mathbf{h_{\{1,2\}}}$ are compact descriptions for the bath Hamiltonians associated with states $\left|1\right\rangle$ and $\left|2\right\rangle$.

There are typically a large number of bath modes present but we are chiefly interested in the diabatic populations and not the minutiae of individual bath coordinates. It is therefore reasonable to treat the bath oscillators as a continuum via a spectral density $J(\omega)$ of the formNitzan (2006):

[TABLE]

It is known that the dynamics of the spin-boson model is completely specified by $V,\epsilon$ and $J(\omega)$Chakravarty and Leggett (1984); Leggett et al. (1987); Nitzan (2006)-though an actual closed form expression remains elusive.

We also note that the spin-boson model is extremely simplified and neglects many physical effects (non-Condon effectsIzmaylov et al. (2011); Mavros, Hait, and Van Voorhis (2016), Duschinski rotationsTang, Lee, and Lin (2003) and bath anharmonicity, to name a few) that are relevant in many condensed phase processes of interest. However, the relative difficulty in solving even this simplified Hamiltonian makes it reasonable to use the spin-boson model as a starting point for method development and then generalize to more complex models.

III Memory Kernels

Memory kernelsNakajima (1958); Zwanzig (1960, 1961, 1964) represent one possible approach for studying the spin-boson and related models. We employ a generalized version of the memory kernel formalism of Sparpaglione and MukamelSparpaglione and Mukamel (1988) in this paper, where the memory kernels $K_{11/22}(t,t_{1})$ control the dynamics of the population $p_{i}(t)$ in the state $\left|i\right\rangle$ in the following manner:

[TABLE]

However, the only independent quantity here is the difference $p_{1}(t)-p_{2}(t)$ as $p_{1}(t)+p_{2}(t)=1$ for all times. This difference is termed as $\left\langle\sigma_{z}\right\rangle(t)$ since $p_{1}(t)-p_{2}(t)=\mathbf{Tr}\left[\bm{\rho}(t)\left(\left|1\right\rangle\left\langle 1\right|-\left|2\right\rangle\left\langle 2\right|\right)\right]=\mathbf{Tr}\left[\bm{\rho}(t)\sigma_{z}\right]$, where $\bm{\rho}(t)$ is the time-dependent density matrix for the whole system. We thus actually just have a single equation:

[TABLE]

where $K_{\pm}(t,t_{1})=K_{11}(t,t_{1})\pm K_{22}(t,t_{1})$.

Time-dependent perturbation theory allows us to expand the population difference $\left\langle\sigma_{z}\right\rangle(t)$ and the kernels $K_{\pm}(t,t_{1})$ into:

[TABLE]

where $\left\langle\sigma_{z}\right\rangle^{(m)}(t)$ and $K^{(m)}_{\pm}(t,t_{1})$ can be found via the procedure described in Appendix A.

The popular non-interacting blip approximation (NIBA) Chakravarty and Leggett (1984); Leggett et al. (1987); Aslangul, Pottier, and Saint-James (1985) approximates $K$ by only using the second-order term (i.e. $K\approx V^{2}K^{(2)}$), yielding quite accurate dynamics in the small $V$ (‘outer sphere’) regime and reducing to Marcus theory for slow bath frequencies at high temperaturesXu and Schulten (1994). It is expected that the dynamics obtained from kernels incorporating higher order terms will be even more accurate, but would be much more computationally expensive to obtain as kernels accurate to $2n$th order in $V$ require evaluation of $2n-1$ dimensional oscillatory integrals. Therefore, alternative approaches that can obtain reasonable approximations of higher order terms from only second-order terms like $\left\langle\sigma_{z}\right\rangle^{(2)}$ are desirable as they would incur a much lower computational cost.

III.1 Initial conditions

The rest of the paper will assume that the initial density matrix $\bm{\rho}(0)=p_{1}(0)\left|1\right\rangle\left\langle 1\right|\otimes\rho_{B}+p_{2}(0)\left|2\right\rangle\left\langle 2\right|\otimes\rho_{B^{\prime}}$, which causes odd order power series terms $\left\langle\sigma_{z}\right\rangle^{(2k+1)}(t)$ and $K_{\pm}^{(2k+1)}(t)$ to be [math]. We will perform relevant derivations with the even simpler initial condition of $\bm{\rho}(0)=\left|1\right\rangle\left\langle 1\right|\otimes\rho_{B}\implies\left\langle\sigma_{z}\right\rangle(0)=1$ for mathematical simplicity, as the expressions derived from this form can be trivially generalized to any initial separable diabatic density matrix due to linearity of time evolution. In particular, we would like to emphasize that Eqn 27 and beyond are true for all $\bm{\rho}(0)=p_{1}(0)\left|1\right\rangle\left\langle 1\right|\otimes\rho_{B}+p_{2}(0)\left|2\right\rangle\left\langle 2\right|\otimes\rho_{B^{\prime}}$.

The initial bath density matrix $\rho_{B}$ is determined by the nature of the dynamical process of interest. Ground state electron transfer processes often use equilibrium $\rho_{B}=\dfrac{e^{-\beta\mathbf{h_{1}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{1}}}\right]}$ ($\beta$ being the inverse temperature $\dfrac{1}{k_{B}T}$) as this describes a scenario where the bath modes are in thermal equilibrium with the reduced Hamiltonian $\mathbf{h_{1}}$ associated with the state $\left|1\right\rangle$ (which has all the population). On the other hand, non-equilibrium $\rho_{B}$ of the type $\dfrac{e^{-\beta\mathbf{h_{2}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{2}}}\right]}$ can be useful for photochemistry, as this corresponds to a situation where the population was purely in $\left|2\right\rangle$ with the bath in thermal equilibrium, before an excitation at $t=0^{-}$ led transfer of all the diabatic population to $\left|1\right\rangle$ without giving the bath modes an opportunity to relax accordingly, due to difference in timescales. We call this initial condition the “photochemical” non-equilibrium condition throughout the rest of the paper for convenience, although it need not correspond to experimentally observable photochemistry for all systems. Other non-equilibrium initial conditions may also make physical sense for different processes. Unless specified otherwise, all our expressions are valid for any $\rho_{B}=\dfrac{e^{-\beta\mathbf{h_{3}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{3}}}\right]}$ where $\mathbf{h_{3}}=\mathbf{h_{0}}+\displaystyle\sum\limits_{i}d_{i}x_{i}$ ($\{d_{i}\}$ can be any arbitrary real number) is a bath Hamiltonian with the same frequencies and modes as $\mathbf{h_{1/2}}$, but with any arbitrary set of equilibrium positions.

III.2 Power series terms

Evaluation of the non-zero even order terms in the power series of $\left\langle\sigma_{z}\right\rangle(t)$ involves tracing over bath modes. Time-dependent perturbation theory reveals that the power series coefficients $\left\langle\sigma_{z}\right\rangle^{(2n)}(t)$ for $\bm{\rho}(0)=\left|1\right\rangle\left\langle 1\right|\otimes\rho_{B}$ are completely specified by integrals of traces $f_{2n}(t_{1},t_{2},t_{3}\ldots t_{2n})$ where:

[TABLE]

For instance, we have:

[TABLE]

Similarly, $\left\langle\sigma_{z}\right\rangle^{(4)}(t)$ is expressed in terms of integrals of $f_{4}$, $\left\langle\sigma_{z}\right\rangle^{(6)}(t)$ by integrals of $f_{6}$ and so on.

III.3 Non-Interacting Blip Approximation

Going up to second order alone, we discover that $K_{11}^{(2)}(t,t_{1})=2\mathbf{Re}\left[f_{2}(t,t_{1})\right]$ is consistent with Eqn. 12. This expression is similar to the non-equilibrium golden rule in Ref [Coalson, Evans, and Nitzan, 1994], although the formalism employed therein is convolution free unlike Eqns 5-6. For the case of $\rho_{B}=\dfrac{e^{-\beta\mathbf{h_{1}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{1}}}\right]}$ (thermal equilibrium initial conditions), this reduces to the well-known Non-Interacting Blip Approximation (NIBA)Caldeira and Leggett (1983); Chakravarty and Leggett (1984); Aslangul, Pottier, and Saint-James (1985). The second-order kernels $K^{(2)}_{11/22}$ are then given byMavros and Van Voorhis (2014):

[TABLE]

where $J(\omega)$ is the spectral density defined earlier in Eqn. 4. It can quickly be noted that these kernels only depend on the time difference $t-t_{1}$ and, as such, we can simply replace the function of two variables $K_{11/22}^{(2)}(t,t_{1})$ with the univariate $K_{11/22}^{(2)}(t-t_{1})$. The resulting kernel $K_{11/22}(t)\approx V^{2}K_{11/22}^{(2)}(t)$ is termed as the NIBA kernel. This approach is expected to yield decent dynamics for small $V$, as it is accurate to the lowest non-zero order for this model.

For non-equilibrium initial conditions however, the traces involved are not time-translationally invariant as they depend on both $(t,t_{1})$ and not just the difference $t-t_{1}$, unlike the equilibrium case. As an example, this generalization for the case of non-equilibrium initial $\rho_{B}=\dfrac{e^{-\beta\mathbf{h_{2}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{2}}}\right]}$results in kernelsMavros, Hait, and Van Voorhis (2016):

[TABLE]

which are very similar to the equilibrium kernels, aside from the subtle difference stemming from the phase $\phi(t,t_{1})$ not being time-translationally invariant. The resulting kernels $K_{11/22}(t,t_{1})\approx V^{2}K_{11/22}^{(2)}(t,t_{1})$ cannot strictly be called NIBA kernels, but this is effectively a natural extension of the idea behind NIBA, and is used in its place for non-equilibrium initial conditions. Similar expressions had been derived earlier for both memory kernelThoss, Wang, and Miller (2001) and non-equilibrium golden rule studiesIzmaylov et al. (2011).

III.4 Resummed Kernels

NIBA is only correct to second order, and performs quite poorly for systems with large diabatic coupling, large energy bias or low temperature. The natural step forward seems to be adding fourth or sixth order terms to the approximate kernel. However, simply stopping there might be problematic as finite truncations of perturbation theory series are often not convergentDyson (1952). Resummation techniques help ameliorate this complication by approximating the unknown higher order terms from the known low order terms. The recipe for memory kernel resummation has been outlined in Refs [Mavros and Van Voorhis, 2014,Gong et al., 2015] and other places, and we will not discuss the details here. We will merely note that this approach yields significantly superior results compared to NIBAMavros and Van Voorhis (2014); Gong et al. (2015); Chen, Berkelbach, and Reichman (2016) but this increased accuracy comes with the much higher computational cost of evaluating $K^{(4)}$ or other higher order kernels exactly.

IV Trace relations

Closed-form expressions for the traces $f_{2n}$ for $n\geq 1$ are necessary for calculating $\left\langle\sigma_{z}\right\rangle^{(2n)}$ or $K^{(2n)}$. We discovered that all traces $f_{2n}$ can be expressed in terms of the second order trace $f_{2}$ alone for any integer $n$ and any choice of the spectral density; if the initial bath density matrix $\rho_{B}$ equals $\dfrac{e^{-\beta\mathbf{h_{3}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{3}}}\right]}$ for $\mathbf{h_{3}}=\mathbf{h_{1}}+\displaystyle\sum\limits_{i}d_{i}x_{i}$ ($\{d_{i}\}$ can be any arbitrary set of real numbers of cardinality $N$). Specifically, we have:

[TABLE]

For example:

[TABLE]

and similarly for $f_{6},f_{8}$ etc.

Essentially, an initial bath density matrix corresponding to a thermal state of any harmonic bath $\mathbf{h_{3}}$ with same frequencies and modes as $\mathbf{h_{1/2}}$ is sufficient to satisfy the trace factorization relationship Eqn. 20 (irrespective of the equilibrium positions of the bath modes in $\mathbf{h_{3}}$) for any bath spectral density. Similar relationships had been presented earlier for equilibrium initial conditions of $\mathbf{h_{3}}=\mathbf{h_{1}}$Skinner and Hsu (1986); Leggett et al. (1987) but we believe that the case of general $\mathbf{h_{3}}$ has not previously been reported. This relationship has been analytically verified to only 12th order for $\mathbf{h_{3}}=\mathbf{h_{1}}$ (equilibrium initial conditions) and 6th order for $\mathbf{h_{3}}=\mathbf{h_{2}}$ (“photochemical” non-equilibrium initial conditions) due to the rather steep memory cost for the evaluation of the integrals needed to find the analytic expressions with MathematicaWolfram (1999). It has however been verified numerically up to $f_{60}$ for multiple $\mathbf{h_{3}}$ resulting from random $\{d_{i}\}$ over a large number of randomly selected time indices $\{t_{i}\}$, indicating that it is correct to 60th order at least, and very likely beyond as well (though as of now, we do not possess a proof for this conjecture). This result is remarkable because it implies that for any spectral density, the dynamics is a function of $f_{2}$ alone. That is to say that, in principle, one can predict all of the dynamics using only $f_{2}(t_{1},t_{2})$ as an input. In analogy with time density functional theory Runge and Gross (1984) (in which the dynamics of a many particle system is written as a functional of the one particle density alone), one can envision a $f_{2}(t_{1},t_{2})$ functional theory - in which the dynamics are directly predicted by some functional of $f_{2}$ alone. NIBA is one such theory, but the result above suggests that the exact result can in principle be constructed from $f_{2}$ alone.

V Case of the Slow bath

Though Eqn. 20 is quite useful in finding all the perturbation theory traces, it does not eliminate the need to numerically integrate $f_{2n}$ over the $2n$ time indices in order to find $\left\langle\sigma_{z}\right\rangle^{(2n)}(t)$. This quickly becomes intractable due to the grid size growing exponentially with $n$, creating a practical upper limit for applying naive perturbation theory even though we can find traces to arbitrary order $2n$ rather simply via Eqn. 20.

However at short times, we may expand $\ln f_{2}(t_{1},t_{2})$ as a Taylor series about $t_{1}=t_{2}=0$ to obtain:

[TABLE]

where

[TABLE]

Higher order terms can be neglected when $\omega t_{1,2}\ll 1$, where $\omega$ is the characteristic response frequency of the system. Taking this short-time/slow-bath approximation and defining a Gaussian function $g(t)=e^{-at^{2}+ibt}$, we obtain:

[TABLE]

which enables us to approximate all these traces as Gaussian functions that can be analytically integrated. Carrying out the integrations over the time indices $t_{i}$ and summing over all orders in $V$ (in the manner described in Appendix C), we find that:

[TABLE]

where $J_{1}$ is a Bessel function of the first kind. We note that Eqn 27 holds for any initial density matrix of the diagonal form $p_{1}(0)\left|1\right\rangle\left\langle 1\right|\otimes\rho_{B}+p_{2}(0)\left|2\right\rangle\left\langle 2\right|\otimes\rho_{B^{\prime}}$, courtesy the linearity of time evolution.

For equilibrium initial conditions we have $b=\epsilon+\lambda$ and $a=\dfrac{2}{\pi}\displaystyle\int\limits_{0}^{\infty}J(\omega)\cosh\dfrac{\beta\omega}{2}\,d\omega$; where $\lambda=\dfrac{4}{\pi}\displaystyle\int\limits_{0}^{\infty}\dfrac{J(\omega)}{\omega}\,d\omega$ is the Marcusian reorganization energy. $a\approx\dfrac{\lambda}{\beta}$ in the limit of frequencies $\{\omega_{i}\}\to 0$. The expression for $\left\langle\sigma_{z}\right\rangle(t)$ obtained with these values for $a,b$ appears to be consistent with the results obtained for a slow bath by Hornbach and DakhnovskiHornbach, Dakhnovskii et al. (1999), which makes us more confident in overall accuracy of our result and the applicability to general initial conditions beyond the equilibrium case. The “photochemical” non-equilibrium initial condition, on the other hand, yields the same $a$, but $b=\epsilon-\lambda$ instead, indicating that for slow-baths at least, equilibrium and “photochemical” conditions lead to same dynamics if $\epsilon=0$.

We can also substitute the exact $\dot{\left\langle\sigma_{z}\right\rangle}^{(2)}(t)=-4\displaystyle\int\limits_{0}^{t}\mathbf{Re}\left[f_{2}(t,t_{1})\right]\,dt_{1}$ in Eqn 27 instead of the Gaussian approximate form given in Eqn 28. This is done in the hope of generalizing Eqn 27 to baths with non-Gaussian $f_{2}$ as the substitution makes Eqn 27 exact to order $V^{2}$ for any set of parameters. We call this version of Eqn 27 GaussSB as it is exact for spin-boson problems with Gaussian $f_{2}$.

VI Benchmarking Accuracy of GaussSB

The lack of a rigorous proof for Eqn. 20 raises some well warranted questions about the accuracy of GaussSB even in the Gaussian limit, leading us to compare it against HEOMIshizaki and Tanimura (2005); Tanimura (2006)—a method known to be numerically exact for $J(\omega)=\dfrac{\lambda}{2}\left(\dfrac{\omega\omega_{c}}{\omega^{2}+\omega_{c}^{2}}\right)$ if a sufficiently large number of Matsubara frequencies are employed and a high hierarchy depth used for calculations. This particular form of $J(\omega)$ appears to be suboptimal for GaussSB on account of the long, slowly decaying high-frequency tail, but we nonetheless stick with it for convenience in benchmarking. We also compared GaussSB to NIBA in order to determine the differences between the two, since they both attempt to approximate the dynamics with the second order term $\left\langle\sigma_{z}\right\rangle^{(2)}(t)$ alone. Comparisons were done in both the slow bath limit (where $f_{2}$ is Gaussian and GaussSB is expected to be exact) and the fast bath limit (where $f_{2}$ is decidedly non-Gaussian and GaussSB is likely to fail).

HEOM calculations were done with the PHI codeStrümpfer and Schulten (2012) while both GaussSB and NIBA were implemented in C++ employing the error function implementation by Johnson Johnson (2012) and the GSL implementation of Quadpack integration routines Gough (2009); Piessens et al. (2012).

VI.1 Slow Baths

GaussSB is expected to be valid in the $\omega_{c}\to 0$ limit (slow bath limit) and so $\omega_{c}$ was set to $10^{-4}$ for all calculations in this subsection. We also set $V=1$ throughout in order to fix the timescale of oscillations in diabatic populations. We first checked for cases without bias ($\epsilon=0$) over a wide range of temperatures and reorganization energies ($0.1\leq\lambda,\beta\leq 10$). Only equilibrium initial conditions were tested, as “photochemical” initial conditions give the same Gaussian parameters $a,b$ for cases where $\epsilon=0$, and thus should yield same dynamics. Four different behaviors were observed in this regime and these are depicted in Fig 1. There is visual agreement between GaussSB and HEOM in each case presented (as well as in many other cases tested that resulted in qualitatively similar dynamics), while NIBA often even fails to reproduce the qualitative behavior. These tests are consistent with our claim that GaussSB is exact in the slow bath limit.

We next consider cases with bias, where equilibrium and “photochemical” dynamics are expected to give different results. Two selected cases are depicted in Fig 2 and we observed similar behavior over other ranges of parameters as well. NIBA is known to be problematic in cases with large $\left|\beta\epsilon\right|$, often yielding absurd $\left|\left\langle\sigma_{z}\right\rangle(t)\right|>1$. We observe the same behavior here, but also find excellent agreement between GaussSB and HEOM, despite GaussSB and NIBA both only employing second order information $\left\langle\sigma_{z}\right\rangle^{(2)}(t)$ alone.

VI.2 Behavior for fast baths

GaussSB cannot be expected to work exactly for fast baths since we explicitly used a small-time approximation taylor series to generate Gaussian traces and sum the perturbation theory series. Numerical tests indicate that this is indeed the case, with GaussSB being qualitatively incorrect in the dissipative regime. We believe that the failure of GaussSB stems from slow decay of coherent oscillations courtesy of the Gaussian approximation (which implicitly assumes slow bath response), as it is the bath response that ultimately causes the transition from oscillatory behavior to dissipation. A related consequence of this is that the formal rate constant associated with GaussSB (obtained from integrating the memory kernel from $t=0$ to $t\to\infty$) is always zero, indicating that GaussSB cannot reproduce correct equilibrium populations or replicate long time exponential decay of the type predicted by Marcus theory and experimentally observed in chemical processes at high temperatures. This represents quite a major obstacle to using pure GaussSB for realistic problems where the dissipative regime is involved.

On the other hand, we do obtain qualitatively correct behavior with GaussSB in regimes where the dynamics is controlled by coherent oscillations as opposed to dissipation (namely cases with large $\beta,V,\epsilon$), but the agreement is not quantitative. We do however note that GaussSB does much better than NIBA in this regime, again despite both of them using only second order information alone. Examples for behavior in both regimes are presented in Fig. 3.

VII Interpolated Kernels

Despite the issues related to applicability to fast baths, we know GaussSB is accurate for time-scales shorter than bath response, as evidenced by its accuracy in the slow bath limit. NIBA is simultaneously known to be accurate at times significantly larger than the bath response as it can reproduce dissipative behavior, indicating that an approach that smoothly ties these two regimes together could potentially be more effective than either. This interpolation scheme should recover pure GaussSB type behavior as $t\to 0$, NIBA type behavior as $t\to\infty$ and behavior intermediate between the two at all times in between.This is distinct from hybrid quantum-classical methods like MC-TDHWang, Thoss, and Miller (2001); Thoss, Wang, and Miller (2001), as it combines exact quantum mechanical results at two limiting cases, as opposed to discriminating between different bath modes by treating only a fraction quantum mechanically.

One issue in constructing an interpolating scheme is that GaussSB and NIBA approach population dynamics in completely different ways. NIBA employs a non-Markovian memory kernel to connect $\dot{\left\langle\sigma_{z}\right\rangle}$ and ${\left\langle\sigma_{z}\right\rangle}$ while GaussSB directly obtains $\dot{\left\langle\sigma_{z}\right\rangle}$ via Eqn 27. The memory kernel typically decays much more rapidly than the population transfer rateKelly et al. (2016), and thus it seemed more appropriate to connect the memory kernels of GaussSB and NIBA together in order to better utilize the short-time exactness of GaussSB. This led us to numerically invert GaussSB ${\left\langle\sigma_{z}\right\rangle}$ to obtain kernels $K_{11/22}^{GSB}(t)$ and then smoothly interpolate between these kernels and the exact second order kernels $K^{(2)}_{11/22}(t,t_{1})$ to get hybrid kernels that tend to pure second-order at long times and pure GaussSB at short times. Mathematically, this implies:

[TABLE]

where $K^{I}$ is the interpolated kernel and $u(t)$ is an interpolating function with the following properties:

$u(t=0)=1$. 2. 2.

$\lim\limits_{t\to\infty}u(t)=0$. 3. 3.

The rate of decay of $u(t)$ (or equivalently, the growth rate of $1-u(t)$) is inversely related to the bath response time. This ensures that GaussSB behavior dies on the timescale of the bath response, since it becomes increasingly inaccurate at long times. 4. 4.

$u(t)$ is well-behaved (i.e., non divergent and continuous, preferably differentiable).

The kernel $K^{I}_{11/22}(t,t_{1})$ is thus almost pure GaussSB when the two time indices $t,t_{1}$ are close, almost pure second order kernel when $t,t_{1}$ are far and intermediate between the two extremes at intermediate separations. For the equilibrium case, $K^{(2)}_{11/22}(t,t_{1})$ is the NIBA kernel, and NIBA behavior is recovered at long times.

There exist many possibilities for $u(t)$, depending on what choices are made regarding bath response times. For the previously employed $J(\omega)=\dfrac{\lambda}{2}\dfrac{\omega\omega_{c}}{\omega^{2}+\omega_{c}^{2}}$ form, it is at least possible to unambiguously define a timescale by means of $\omega_{c}$, but such an option is not open for all spectral densities. This leads us to obtain decay rate via $\theta(t)=\dfrac{\displaystyle\int\limits_{0}^{\infty}J(\omega)\cos\omega t\coth\dfrac{\beta\omega}{2}\,d\omega}{\displaystyle\int\limits_{0}^{\infty}J(\omega)\coth\dfrac{\beta\omega}{2}\,d\omega}$ as $\theta(t)$ corresponds to the normalized energy gap fluctuation autocorrelation function $\dfrac{\mathbf{Re}\left[\left\langle\delta\Delta E(t)\delta\Delta E(0)\right\rangle\right]}{\left\langle\delta\Delta E(0)\delta\Delta E(0)\right\rangle}$, where the numerator is a quantity often computed in MD simulations to obtain $J(\omega)$ and is thus readily available in most cases. This decays smoothly from $1\to 0$ as $t$ runs from $0\to\infty$. Consequently, $\theta(t)$ is a measure of the dissipation rate of initial energy fluctuations and thus gives a timescale for switching to the dissipative regime. Furthermore, $\theta(t)=e^{-\omega_{c}t}$ for Debye spectral densities at small $\beta$ (high temperatures, i.e. the classical limit), indicating that it in fact decays with the same timescale as the bath frequencies.

$\theta(t)$ thus is a good metric for the rate of decay. That is however not the only factor we should consider, as although $\theta(t)$ contains information about bath frequencies as well as some temperature effects, it does not at all account for the energy gap $\epsilon$ which is known to affect population coherence time-scales as well. Furthermore, NIBA is known to be especially terrible for large $\left|\beta\epsilon\right|$ (as multiple figures in the preceding section highlight) and this weakness also ought to be considered for building $u(t)$. We therefore propose $u(t)=\left|\theta(t)\right|^{p}$ to be a general form for the interpolating function $u(t)$, where $p$ is very small for $\beta\epsilon\gg 1$ (allowing GaussSB to dominate) but reduces to just $1$ for $\epsilon=0$. $p=\mathrm{sech}\beta\epsilon$ was chosen by us as it satsifies the above requirements and corresponds to the $e^{\beta A\left(\epsilon\right)}$, where $A\left(\epsilon\right)=-\dfrac{1}{\beta}\ln\cosh\beta\epsilon$ is the component of the Gibbs-Boguliubov free energy dealing with biasSong, Wang, and Van Voorhis (2008). NIBA performance is known to deteriorate as $A(\epsilon)$ increases in magnitudeSong, Wang, and Van Voorhis (2008), indicating that $p$ can be a decent metric for controlling the GaussSB to NIBA switching rate. The absolute value in $u(t)=\left|\theta(t)\right|^{p}$ is taken since though $\theta(t)\geq 0$ in most cases, a few pathological parameter sets might lead to $\theta(t)<0$, which would make taking non-integer powers problematic. Taking $\left|\theta(t)\right|$ solves this issue, though it can lead to non-differentiable $u(t)$ at points where $\theta(t)=0$ changes sign. This is unfortunate, but we feel it is an acceptable compromise due to its rarity and because it does not lead to non-differentiable $\left\langle\sigma_{z}\right\rangle$ as it only affects the kernel.

For determining the accuracy of this interpolated approach, we considered a few Hamiltonian parameters studied earlier via memory kernel resummation in Ref [Mavros and Van Voorhis, 2014] (all using equilibrium initial conditions) and compared the performance of HEOM, NIBA, Interpolated GaussSB and Padé resummed memory kernelsBaker and Graves-Morris (1996) (obtained from data generated for Ref [Mavros and Van Voorhis, 2014]) correct to fourth order in $V$. The results for cases with $\epsilon=0$ are given in Fig 4 and it appears that the interpolated kernel approach is comparable to (and for large $V$ is superior to) dynamics obtained from the much more computationally expensive fourth-order resummed kernel approach. It thus appears that the interpolated kernel approach could be a cheap and effective way to incorporate $V^{4}$ and higher order terms to NIBA at a significantly lower computational cost than numerical evaluation of higher order kernels for resummation. We also note that in addition to our chosen form of $u(t)$, we experimented with some different forms of $u(t)$ (decaying exponentials and Gaussians in $\omega_{c}t$) for these cases and obtained roughly similar behavior. This aspect however was not explored in greater detail as other $u(t)$ explored did not easily generalize to non Drude-Lorentz spectral densities where a characteristic $\omega_{c}$ is not apparent.

Cases with a non-zero $\epsilon$ pose a greater challenge, as we know that neither GaussSB nor NIBA can accurately determine the equilibrium populations in such cases and thus the combined approach is unlikely to succeed in this task either. It has been shownMavros and Van Voorhis (2014) that it is possible to adjust the resummation scheme to yield kernels correct to fourth order that recover equilibrium populations in a least-squares sense. This indicates that the interpolated kernel in these cases should be expected to perform worse than the optimal fourth order resummed kernel as the latter has more information than the former. Fig 5 depicts the comparison of HEOM, NIBA, intepolated kernel and optimized Landau–Zener resummationSumi (1980) (found to be superior to Padé for cases with biasMavros and Van Voorhis (2014)) for some systems studied in Ref [Mavros and Van Voorhis, 2014]. This comparison indicates that the interpolated kernel approach gives dynamics of comparable quality to fourth-order resummed kernels, although it tends to perform worse at long times as it does not contain information about $\left\langle\sigma_{z}\right\rangle(t\to\infty)$ that the resummation technique employs. Nonetheless, it recovers quite reasonable short-time dynamics, with less information and at a fraction of the cost relative to the fourth order methods. Performance for the difficult cases where exact $\left\langle\sigma_{z}\right\rangle(t\to\infty)$ is close to $\pm 1$ (such as Fig 5(a)) can potentially be further improved by altering $u(t)$ to ensure recovery of correct equilibrium populations (perhaps via tuning the power parameter $p$).

VIII Conclusion and Future Directions

We have found a relationship between perturbation theory traces of the spin-boson model that allowed us to solve for population dynamics analytically for the case of a Gaussian trace $f_{2}(t_{1},t_{2})$. We call this method GaussSB and validate its accuracy by benchmarking it against the highly accurate HEOM method for infinitely slow baths. GaussSB is unfortunately not particularly useful for many problems of chemical interest since $f_{2}(t_{1},t_{2})$ is only Gaussian for baths with small frequencies, though it appears to recover the non-dissipative behavior at low temperatures with qualitative accuracy for even fast baths.

We remedy this by formulating a hybrid kernel that is exact to second order in $V$ and includes higher order correction terms incorporated via GaussSB that decay on the timescale of the bath. This decay is controlled by an interpolating function $u(t)$, in order to prevent higher order GaussSB terms (which become increasingly non-exact at long times) from becoming too influential on timescales larger than bath response. This hybrid approach yields qualitatively accurate dynamics for even fast baths that is comparable in quality to results obtained from using more computationally expensive resummed memory kernels accurate to fourth order. The hybrid kernel method is especially effective in regimes where NIBA normally breaks down (large $\beta,V,\epsilon$), indicating that the interpolation approach was successful in incorporating critical amounts of higher order behavior.

We have thus devised a hybrid method that interpolates between NIBA and GaussSB to give quite accurate dynamics for the spin-boson model for the same computational complexity as NIBA (since it does not need evaluation of expensive triple or higher order integrals). This hybrid approach could thus be used in place of NIBA to explore dynamics across a large parameter space since it has superior performance to NIBA in traditionally problematic regimes and is not observed to ever perform any worse that NIBA, on account of it incorporating higher order effects that NIBA misses.

It should in principle be possible to generalize this approach to other model Hamiltonians like the linear vibronic coupling (LVC) modelMavros, Hait, and Van Voorhis (2016), as long as equivalent slow-bath solutions like GaussSB exist. In practice, there are problems stemming from the fact that simple trace factorization relations like Eqn 20 have not been found (and may simply not exist) for more sophisticated systems. Directly employing the GaussSB kernel (using the exact $\left\langle\dot{\sigma_{z}}\right\rangle^{(2)}(t)$ for the model) might be an acceptable compromise for such cases, even though it would sacrifice exactness in the slow-bath limit. The resulting method would resemble a resummation scheme as it would employ the exact second order term to approximate higher order terms, and would furthermore be much cheaper than actually computing fourth order kernels for the problem. We hope to explore the effectiveness of this approach in the future to determine whether it has any promise.

In the future, we also seek to explore routes to generalize GaussSB to have kernels correct to fourth (or higher) order at all times and see if that yields improved dynamics. The easiest way to achieve this is via a generalization of Eqn 27 of the form

[TABLE]

where the function $f(t_{1})$ is chosen to ensure $\dot{\left\langle\sigma_{z}\right\rangle}(t)$ is exact to fourth power in $V$. The kernel resulting from this expression will be exact to fourth order as it contains both second and fourth order information. Generalizing even further, it appears that it might possible to generate other resummation methods employing the general form of Eqn. 27 that are exact in the limit of Gaussian $f_{2}$. They would still incur significant computational expense in order to have exact fourth (or higher) order behavior, but could potentially prove more effective than pure Padé or Landau–Zener resummation and might not even need interpolation on account of the extra information supplied. It may also be possible to recover the correct equilibrium populations with such kernels, and we hope to study the accuracy of such kernels (both with interpolation and without) in the future to determine whether this is indeed the case.

Acknowledgements.

This work was supported by NSF Grant No. CHE-1058219. D.H. would like to thank the MIT UROP office for support. M.G.M. would also like to thank the NSF GRFP program for funding.

Appendix A Time-Dependent Perturbation Theory

Time-dependent perturbation theory is the necessary first step for any perturbative approaches to quantum dynamics. We treat the diabatic coupling $V$ as the perturbation, which allows us to define the zeroth order Hamiltonian , consequently enabling us to express operator $A(t)$ in the interaction picture as $A_{I}(t)=e^{i\mathbf{H_{0}}t}A(t)e^{-i\mathbf{H_{0}}t}$. The perturbation in the interaction picture is thus $V_{I}(t)=e^{i\mathbf{H_{0}}t}\begin{pmatrix}0&V\\ V&0\end{pmatrix}e^{-i\mathbf{H_{0}}t}=V\sigma^{x}_{I}(t)$ ($\sigma^{x}$ being the Pauli matrix in the x direction) and the interaction picture density matrix can be expressed as $\bm{\rho}_{I}(t)=e^{i\mathbf{H_{0}}t}\bm{\rho}(t)e^{-i\mathbf{H_{0}}t}$ (where $\rho(t)$ is the actual time dependent density matrix). The von-Neumann equation for the time-evolution of the interaction-picture density matrix allows us to state that:

[TABLE]

where $\bm{\rho}^{(m)}_{I}(t)$ are independent of $V$.

The diabatic populations $p_{1/2}(t)$ can be obtained from the density matrix by tracing out the bath modes and taking the diagonal elements. Mathematically:

[TABLE]

$p_{1/2}(t)$ are not independent as $p_{1}(t)+p_{2}(t)=1$, and thus the population dynamics can be specified with only knowledge of the difference $p_{1}(t)-p_{2}(t)$. This is termed as $\left\langle\sigma_{z}\right\rangle(t)$ as $p_{1}(t)-p_{2}(t)=\mathbf{Tr}\left[\bm{\rho}_{I}(t)\left(\left|1\right\rangle\left\langle 1\right|-\left|2\right\rangle\left\langle 2\right|\right)\right]=\mathbf{Tr}\left[\bm{\rho}_{I}(t)\sigma_{z}\right]$.

Thus we can express $\left\langle\sigma_{z}\right\rangle(t)$ as a power series in $V$, with the $m$th order term $\left\langle\sigma_{z}\right\rangle^{(m)}(t)=\mathbf{Tr}_{B}\left[\bm{\rho}^{(m)}_{I}(t)\sigma_{z}\right]$ for all non-negative integers $m$. Expressing the kernels $K_{11/22}(t)=\displaystyle\sum\limits_{m=0}^{\infty}V^{m}K^{(m)}_{11/22}(t)$ as a power series in $V$ and equating terms with same order of $V$ on both sides of Eqns 5 and 6 allows us to express power series coefficients $K^{(m)}_{11/22}(t)$ in terms of $\left\langle\sigma_{z}\right\rangle^{(m)}(t)$, which we have already shown can be found from perturbation theory. Going up to some finite order $m$ and finding all $\left\langle\sigma_{z}\right\rangle^{(k)}(t)$ for $k\leq m$ would thus provide us with the correct kernel coefficients $K^{(k)}_{11/22}(t)$, which can be directly employed or be resummed to approximate the exact kernel correct to all orders in $V$.

Appendix B Taylor series for $\ln f_{2}(t_{1},t_{2})$

We have:

[TABLE]

Since

[TABLE]

Expanding $\ln f_{2}(t_{1},t_{2})$ as a Taylor series about $t_{1}=t_{2}=0$ gives:

[TABLE]

Higher order terms can be neglected if:

[TABLE]

etc are small, which is the case if the bath frequencies $\{\omega_{i}\}$ are small. As an example, the combined third order term for equilibrium $\rho_{B}=\dfrac{e^{-\beta\mathbf{h_{1}}}}{\mathbf{Tr}\left[e^{-\beta\mathbf{h_{1}}}\right]}$ is $-\dfrac{2i}{3\pi}(t_{1}-t_{2})^{3}\displaystyle\int\limits_{0}^{\infty}\omega J(\omega)\,d\omega$, which is guaranteed to be small if $J(\omega)$ is only considerable for small $\omega$. Similarly, the fourth order term for this case is: $\dfrac{(t_{1}-t_{2})^{4}}{6\pi}\displaystyle\int\limits_{0}^{\infty}\omega^{2}J(\omega)\coth\dfrac{\beta\omega}{2}\,d\omega$, which will be small if $J(\omega)$ is only large for small $\omega$.

Appendix C Derivation of Slow Bath Solution

Integrating Gaussian $f_{2n}$ analytically yields a relation between population growth terms $\dot{\left\langle\sigma_{z}\right\rangle}^{(2n)}$:

[TABLE]

for all integers $n>1$. This form is more convenient to handle in Laplace space where we have:

[TABLE]

by repeated application of Eqn 51. This permits us to sum the infinite power series for $\dot{\left\langle\sigma_{z}\right\rangle}$ in Laplace space, as we have:

[TABLE]

Returning to the time domain, we have via known inverse Laplace transformsErdélyi and Bateman (1954):

[TABLE]

where $a\neq 0$. $J_{1}$ is a Bessel function of the first kind, giving us an analytical solution to the spin-boson model for all baths with trace $f_{2}$ Gaussian in time.

Appendix D GaussSB in the limit of Rabi Oscillations

If we consider the trivial case of $J(\omega)=0$, we obtain $f_{2}(t_{1},t_{2})=e^{i\epsilon(t_{1}-t_{2})}$, which is a complex exponential and not actually a Gaussian. Thus $a=0$ and we can’t use Eqn. 28 directly, but it is easy to use Eqn. 59 and find:

[TABLE]

[TABLE]

which is the exact result that can be directly calculated for this trivial example, showing that it correctly recovers the Rabi limit.

Appendix E Kernels from Populations

GaussSB gives $\left\langle\sigma_{z}\right\rangle(t)$ directly and does not employ any memory kernels. Nonetheless, it would be preferable to interpolate between time non-local memory kernels that population growth rates to fully exploit the the complete short-time information contained in GaussSB, compelling us to determine ways to invert the GaussSB populations $\left\langle\sigma_{z}\right\rangle(t)$ to yield kernels.

We begin from the one time-index version of the Sprapaglione-Mukamel formalismSparpaglione and Mukamel (1988) as there is no way (or reason) to introduce a second time index whose behavior we cannot solely recover from the univariate $\left\langle\sigma_{z}\right\rangle(t)$. Specifically, we use a variant of Eqn 7 where:

[TABLE]

Differentiating Eqn. 66 with respect to time, we obtain:

[TABLE]

Let us consider two initial conditions: $\left\langle\sigma_{z}\right\rangle(0)=\pm 1$ and call the corresponding $\left\langle\sigma_{z}\right\rangle(t)$ as $\left\langle\sigma_{z}\right\rangle_{\pm}(t)$. Then we have:

[TABLE]

GaussSB gives us $\dot{\left\langle\sigma_{z}\right\rangle}_{\pm}(t-t_{1})$ from which we can obtain the second derivatives $\ddot{\left\langle\sigma_{z}\right\rangle}_{\pm}(t-t_{1})$ via direct differentiation of Eqn 27 or finite differences (which was the route we took). Then we can solve Eqn. 69 numerically on a grid for $K_{+}(t)$ (by treating the integral as a finite sum via trapezoid rule and solving the resulting system of linear equations for instance, as we have done: but other approaches like Archimedes summation or other Newton-Coates approaches should be equally valid), making use of the fact that $K_{+}(0)=4V^{2}$ since only second order terms matter then. Finally, back substitution of $K_{+}(t)$ in Eqn 68 is sufficient to recover $K_{-}(t)$, yielding the kernels.

Appendix F Kernels in Laplace Space

Eqn. 69 offers an interesting alternative approach for accessing kernels via Laplace space. Taking Laplace transforms on both sides, we obtain:

[TABLE]

We can determine $\mathcal{L}\left[\dot{\left\langle\sigma_{z}\right\rangle}^{(2)}_{\pm}\right](s)$ via simple Laplace transform from time-domain and then use it to find $\mathcal{L}\left[K_{+}\right](s)$ (and similarly $\mathcal{L}\left[K_{-}\right](s)$). However, inverse Laplace transforms are numerically unstable and thus using this route to obtain time domain kernels is suboptimal.