A Unified Stochastic Formulation of Dissipative Quantum Dynamics. II. Beyond Linear Response of Spin Baths

TL;DR

This paper extends the generalized hierarchical equation of motion to include higher order effects in spin bath models, revealing when linear response approximations are valid or break down due to non-Gaussian and non-Markovian dynamics.

Contribution

It introduces a systematic method to incorporate anharmonic effects in spin bath decoherence models, analyzing the validity of linear response approximations in different physical scenarios.

Findings

Linear response holds for discretized baths with large N_B due to suppression of non-Gaussian effects.

Physical spin baths often violate linear response assumptions due to direct spin-spin interactions.

Symmetries in spin Hamiltonians can lead to non-Markovian dynamics beyond linear response.

Abstract

We use the "generalized hierarchical equation of motion" proposed in Paper I to study decoherence in a system coupled to a spin bath. The present methodology allows a systematic incorporation of higher order anharmonic effects of the bath in dynamical calculations. We investigate the leading order corrections to the linear response approximations for spin bath models. Two types of spin-based environments are considered: (1) a bath of spins discretized from a continuous spectral density and (2) a bath of physical spins such as nuclear or electron spins. The main difference resides with how the bath frequency and the system-bath coupling parameters are chosen to represent an environment. When discretized from a continuous spectral density, the system-bath coupling typically scales as $\sim 1/\sqrt{N_B}$ where $N_B$ is the number of bath spins. This scaling suppresses the non-Gaussian characteristics of the spin bath and justify the linear response approximations in the thermodynamic limit. For the physical spin bath models, system-bath couplings are directly deduced from spin-spin interactions with no reason to obey the $1/\sqrt{N_B}$ scaling. It is not always possible to justify the linear response approximations. Furthermore, if the spin-spin Hamiltonian and/or the bath parameters are highly symmetrical, these additional constraints generate non-Markovian and persistent dynamics that is beyond the linear response treatments.

Figures

Equations

$$\hat{H}$$

$$\displaystyle\hat{H}_{\text{int}}=\left\{\begin{array}[]{l}\sum_{k>0}g_{k}\hat{\sigma}^{z}_{0}\hat{\sigma}^{x}_{k},\\ \sum_{k>0}g_{k}\hat{\sigma}^{z}_{0}\hat{\sigma}^{z}_{k},\\ \sum_{k>0}g_{k}(\hat{\sigma}^{+}_{0}\hat{\sigma}^{-}_{k}+\hat{\sigma}^{-}_{0}\hat{\sigma}^{+}_{k}),\\ \sum_{k>0}g_{k}(\hat{\sigma}^{x}_{0}\hat{\sigma}^{x}_{k}+\hat{\sigma}^{y}_{0}\hat{\sigma}^{y}_{k}+\hat{\sigma}^{z}_{0}\hat{\sigma}^{z}_{k}).\end{array}\right.$$

$${\hat{H}}_B$$

$${\hat{H}}_{\text{int}}$$

$$E_{n,k}={\omega }_k\left(n+\frac{1}{2}\right)-{\chi }_k{\left(n+\frac{1}{2}\right)}^2,$$

$$D_k=\frac{{\omega }_k\Lambda }{2},{\alpha }_k=\sqrt{\frac{{\omega }_k}{\Lambda }},\text{and}{\chi }_k=\frac{{\omega }_k}{2\Lambda },$$

$${\hat{H}}_B$$

$${\hat{H}}_{\text{int}}$$

$${\hat{H}}_B$$

$${\hat{H}}_{\text{int}}$$

$$J_{\text{eff}}(\omega )=\tanh \left(\frac{\beta \omega }{2}\right)J(\omega ).$$

$$d{\tilde{\rho }}_s=-idt\left[{\hat{H}}_s,{\tilde{\rho }}_s\right]-idt\mathcal{B}(t)\left[A,{\tilde{\rho }}_s\right]$$

$$-\frac{i}{\sqrt{2}}d{W}^{\ast }A{\tilde{\rho }}_s+\frac{i}{\sqrt{2}}d{V}^{\ast }{\tilde{\rho }}_{s}A,$$

$$d{\tilde{\rho }}_B=-idt\left[{\hat{H}}_B,{\tilde{\rho }}_B\right]+\frac{1}{\sqrt{2}}dW\left(B-\mathcal{B}(t)\right){\tilde{\rho }}_B$$

$$+\frac{1}{\sqrt{2}}dV{\tilde{\rho }}_B\left(B-\mathcal{B}(t)\right).$$

$$\mathcal{B}(t)=\sum _k{g}_k{\text{Tr}}_B\left\{{\tilde{\rho }}_B(t){\sigma }_k^x\right\}.$$

$${\rho }_s(t)=\mathcal{E}({\tilde{\rho }}_s(t)).$$

$$\mathcal{B}(t)=\frac{1}{\sqrt{2}}\left({\int }_0^{t}dW(s)\alpha (s)+{\int }_0^{t}dV(s){\alpha }^{\ast }(s)\right),$$

$$\mathcal{B}(t)$$

$$\frac{d{\rho }_s}{dt}=-i\left[{\hat{H}}_s,{\rho }_s\right]-i\left[A,\mathcal{E}\left(\mathcal{B}(t){\tilde{\rho }}_s\right)\right].$$

$${\Phi }_{n+1,m}(t,t_1,\dots ,t_n)=\sum _{\mathbf{j}}{\chi }_{\mathbf{j}}^{n+1,m}{\varphi }_{j_1}(t-t_1)\cdots {\varphi }_{j_n}(t_n-t_1),$$

$$\frac{d}{dt}{\varphi }_j(t)=\sum _{j^{\prime }}{\eta }_{j{j}^{\prime }}{\varphi }_{j^{\prime }}(t).$$

$$\displaystyle\mathbf{A}_{n}=\left[\begin{array}[]{l}a^{n}_{1\mathbf{j}_{1}}\cdots a^{n}_{1\mathbf{j}_{k}},0,\dots\\ \vdots\\ a^{n}_{2^{n}\mathbf{k}_{1}}\cdots a^{n}_{2^{n}\mathbf{j}_{k}},0,\dots\end{array}\right],$$

$$a_{m\mathbf{j}}^n\equiv {\left(\frac{1}{\sqrt{2}}\right)}^n{\int }_0^t{\int }_0^{s_1}\cdots {\int }_0^{s_{n-1}}dU(s_1)\dots dU(s_n){\varphi }_{j_1}(t-s_1){\varphi }_{j_2}(s_2-s_1)\cdots {\varphi }_{j_n}(s_n-s_1),$$

$$\mathcal{B}(t)=\sum _{n,m,\mathbf{j}}{\chi }_{\mathbf{j}}^{n+1,m}{a}_{m\mathbf{j}}^n.$$

$${\rho }^{\left[{\mathbf{A}}_1\right]\left[{\mathbf{A}}_2\right]\left[{\mathbf{A}}_3\right]}\cdots \equiv \mathcal{E}\left(\prod _{n,m,k}{a}_{m{\mathbf{j}}_k}^n{\tilde{\rho }}_s(t)\right),$$

$$\mathcal{E}\left(\mathcal{B}(t){\tilde{\rho }}_s\right)=\sum _{n,m,\mathbf{j}}{\chi }_{\mathbf{j}}^{n+1,m}{\rho }^{\dots \left[{\mathbf{A}}_n\right]\dots }(t),$$

$${\partial }_t{\rho }^{\left[{\mathbf{A}}_1\right]\left[{\mathbf{A}}_2\right]\left[{\mathbf{A}}_3\right]}=-i\left[H_s,{\rho }^{\left[{\mathbf{A}}_1\right]\left[{\mathbf{A}}_2\right]\left[{\mathbf{A}}_3\right]}\right]-i\sum _{n,m,\mathbf{j}}{\chi }_{\mathbf{j}}^{n,m}\left[A,{\rho }^{\cdots \left[{\mathbf{A}}_n+(m,\mathbf{j})\right]\cdots }\right]$$

$$-i\sum _{n,m,\mathbf{j}}{\varphi }_{j_1}(0)A{\rho }^{\cdots \left[{\mathbf{A}}_{n-1}+(m^{\prime },{\mathbf{j}}_1)\right]\left[{\mathbf{A}}_n-(m,\mathbf{j})\right]\cdots }-i\sum _{n,m,\mathbf{j}}{\varphi }_{j_1}(0){\rho }^{\cdots \left[{\mathbf{A}}_{n-1}+(m^{\prime },{\mathbf{j}}_1)\right]\left[{\mathbf{A}}_n-(m,\mathbf{j})\right]\cdots }A$$

$$+\sum _{n,m,\mathbf{j}}{\eta }_{\mathbf{j}{\mathbf{j}}^{\prime }}{\rho }^{\cdots \left[a_{m\mathbf{j}}^n\to a_{m{\mathbf{j}}^{\prime }}^n\right]\cdots }$$

$$\hat{\rho }(0)$$

$$⟨\uparrow |{\rho }_s(t)|\downarrow ⟩=⟨\uparrow |{\rho }_s(0)|\downarrow ⟩e^{-iϵt}{e}^{\Gamma (t)},$$

$$\Gamma (t)$$

$${\hat{\sigma }}_k^{\pm }(t)$$

Full text

I Introduction

Understanding the dissipative quantum dynamics of a system embedded in a complex environment is an important topic across various sub-disciplines of physics and chemistry. Significant progress in the understanding of condensed phase dynamics has been achieved within the context of a few prototypical modelsLeggett et al. (1987); Breuer and Petruccione (2002); U (2012) such as Caldeira-Leggett model and spin-boson model. The environment is often modeled as a bosonic bath characterized with a spectral density, from which bath-induced decoherence can be deduced. This prevalent adoption of bosonic baths is based on following considerations: (1) the simple mathematical tractability of a Gaussian bath model, (2) the linear response of an environment is often sufficient to account for quantum dissipations and (3) the spectral density can be extracted from the classical molecular dynamics simulations

Despite the above-mentioned merits, there exists scenarios in which the “bosonization” of an environment is inadequate. For instance, the electron-transfer reaction in condensed phase is often approximated with the spin boson model. The abstract model treats the generic quantum environment as a set of harmonic oscillators, which corresponds to taking only linear response of solvent effects outside a solvation shell while imposing a “harmonic approximation” on the vibrations modes inside the shell. The anharmonicity and higher-order nonlinear response can be substantial when the donor-acceptor complex is strongly coupled to some low-frequency vibrational modes or present in a nonpolar liquids. To better understand the anharmonic effects of the environment, several groups including oursKryvohuz and Cao (2005); Wu and Cao (2001) have studied correlation functions of anharmonic oscillators and a generalized spin boson model with a bath of independent Morse or quartic oscillators. Similarly, a spin bath can be viewed as an extreme limit of anharmonic oscillators and provides additional insights into condensed phase dynamics.

On the other hand, physical spin bath models, corresponding to localized electron and/or nuclear spins, have received increased attentions due to ongoing interests in developing various solid-state quantum technologies Hsieh et al. (2012); Kloeffel and Loss (2013); Gao et al. (2015) under the ultralow temperature regime when interactions with the phonons or vibrational modes are strongly suppressed. For these spin-based environments, the spectral density is no longer a convenient characterization of the bath. Instead, each bath spin is explicitly specified with parameters $\{\omega_{k},g_{k}\}$, the intrinsic energy scale and the system-bath coupling coefficients, respectively.

In this work, we investigate corrections to the standard linear response treatment of quantum dissipations induced by a spin bathProkofiev and Stamp, P. C. E. (2000). To quantitatively capture these higher order responses, we utilize our recently proposed generalized hierarchical equation of motion (gHEOM) methodHsieh and Cao (2016) to incorporate higher order cumulants of an influence functional into an extended HEOM frameworkTang et al. (2016) through a stochastic formulationStockburger and Grabert (2002); Shao (2004); Lacroix (2005). Even though it is possible to derive the gHEOM directly through the path integral influence functional approachTanimura (2006); Makri (1999), we emphasize that the stochastic approach provides an extension to additional methodology developments such as hybrid deterministic-stochastic algorithmsMoix and Cao (2013). This is a direction we are pursuing. Due to the enormous numerical efforts needed to generate accurate long-time results, we restrict the numerical illustrations in the short-time limit. Recently, our groupCerrillo and Cao (2014) and others have proposed methods to construct the memory kernel of a generic bath from numerically exact short-time dynamical results. Hence, the gHEOM provides an invaluable tool to capture the anharmonicity and non-Gaussian effects of a generic quantum environment when used along with these other methods to correctly reproduce long-time results. Furthermore, starting from the stochastic formalism or the gHEOM, it also serves as a starting point to derive master equationsJang et al. (2002) incorporating these higher order nonlinear effects and can be more efficiently solved to extract long-time results.

As alluded earlier, we cover both scenarios: the spin bath as a specific realization of an anharmonic condensed-phase environment and the physical spin bath. In this work, we always explicitly take a spin bath as a collection of finite number of spins. For physical spin bath models, this restriction reflects the reality that there could only be a finite number of spins in the surrounding environment. For an anharmonic condensed phase environment, if we simply take the spin bath as a realization of an infinitely large heat bath then it has been rigorously shownSuárez and Silbey (1991); Makri (1999) that all higher order response function must vanish in the thermodynamic limit. On the other hand, if we perform atomic simulation of solvents in a condensed phase, the anharmonicity can probably be attributed to a few prominent degrees of freedom. Therefore, we restrict the number of bath spins in order to probe the effects of higher order response functions. Many earlier studiesCaldeira et al. (1993); Segal (2014); Wang and Shao (2012); Gelman et al. (2004); Zhang et al. (2015); Lü and Zheng (2009) on the spin bath focused on the thermodynamic limit and restricted to analyze the linear response only. Some interesting phenomena include more coherent population dynamicsWang and Shao (2012) in the nonadiabatic regime at elevated temperature and the onset of negative thermal conductivitySegal (2014) in a molecular junction coupled to two large spin baths held at different temperatures. In App. B, we report our own investigation on differences between a spin bath and a bosonic bath in the linear response limit when the spin bath can be accurately mapped onto an effective bosonic bath.

In this work, the main difference between the two types of environment comes down to how the parameters $\{\omega_{k},g_{k}\}$ are assigned to each bath spin as explained later. In general, we find the anharmonicity is more pronounced at the low-frequency end when the spectral density for condensed phase environment is the commonly used Ohmic form. Therefore, a slow spin bath at low temperature could potentially pose the most difficulty for a linear response treatment of the bath. On the other hand, for a physical spin bath model, highly non-Markovian and persistent dynamicsBortz and Stolze (2007); Wang et al. (2013); Chen et al. (2007); Seifert et al. (2016); Breuer et al. (2004) emerge under a combination of narrowly distributed bath parameters and highly symmetrical system-bath Hamiltonians. To accurately reproduce these results might require taking higher order response of the spin bath into account.

The paper is organized as follows. In Sec.II, we introduce the spin bath models of interest. In Sec.III, we provide a brief account of the stochastic formalismHsieh and Cao (2016) and how to use it to derive the gHEOM with a systematic inclusion of higher order cumulants of an influence functional. In Sec.IV, we first discuss an exactly solvable dephasing model to stress the importance of higher order cumulant corrections and present a benchmark to validate our numerical method then move on to study both finite size representation of the condensed-phase environment and an Ising spin bath. A brief summary is given in Sec.V. In App.A, we provide additional materials on the stochastic derivation of the generalized HEOM. In App.B, we investigate the linear response effects of the spin bath and physical signatures that one can use to distinguish a spin based condensed-phase environment from a bosonic one.

II Spin Bath Models

We consider the following Hamiltonian in this work,

[TABLE]

where $\hat{H}_{B}=\sum_{k>0}(\omega_{k}/2)\hat{\sigma}^{z}_{k}$ and the standard partition of the system (subscript s), bath (subscript B) and the mutual interaction (subscript int) is implied. The general spin-spin interaction takes the form $\hat{\sigma}^{\alpha}_{0}\hat{\sigma}^{\beta}_{1}$, where $\alpha/\beta$ denotes the Cartesian components of Pauli matrices. Among the choices, most common system-bath interactions read,

[TABLE]

In this work, we should explicitly consider first two interaction forms. The first form is appropriate for modelling condensed-phase environment, while the second form is useful in the decoherence study in quantum computing and related contexts.

II.1 Anharmonic condensed-phase environment

One simple-model approach to investigate anharmonicity of a condensed-phase environment is to generalize the typical bosonic bath by substituting harmonic oscillators with the quartic oscillators or Morse oscillatorsLópez-López et al. (2011). We briefly illustrate the steps to obtain a spin bath model corresponding to the low-energy spaces of a bath of Morse oscillatorsKryvohuz and Cao (2005).

We still begin with Eq. (1) but having a different bath part,

[TABLE]

The $n$-th eigen-energy of the $k$-th Morse oscillator is given by

[TABLE]

where the fundamental frequency $\omega_{k}=\alpha_{k}\sqrt{2D_{k}}$ and the anharmonicity factor $\chi_{k}=\alpha^{2}_{k}/2$. Following Ref. López-López et al., 2011, one can characterize the anharmonicity of each mode by imposing a parameter, $\Lambda$ (the number of bound states in a Morse oscillator). Under the condition of fixed $\Lambda$, one gets

[TABLE]

for each Morse oscillator with a free parameter $\omega_{k}$. As clearly implied in Eqs. (II.1) and (8), the Morse potential and the energy level spacing smoothly reduce back to those of a harmonic oscillator in the limit of $\Lambda\rightarrow\infty$. The recovered harmonic bath is characterized by

[TABLE]

On the other hand, by setting $\Lambda=2$, an effective spin bath emerges. The Equation (II.1) can now be cast as

[TABLE]

which correspond to the first interaction form in Eq. (6).

The mapping of a generic anharmonic environment onto a spin bath is more universal than the specific example of Morse oscillators. A general approach to achieve the mapping is to formulate an influence functional of the bath then perform a cumulant expansion, which characterizes the bath-induced decoherence through multi-time correlation functions. One then has a clear set of criteria to construct a spin bath to reproduce the bath’s response up to a specific cumulant expansion order. This is achievable as a set of spins (or qubits) constitute a versatile quantum simulatorLloyd (1996) that can simulate other simple quantum systems.

II.2 Physical Spin Bath

In the present context, the spin bath is not just a conceptual model but represent the actual spin-based environment composed of nuclear / electron spins in the surrounding medium of a physical system. Depending on details regarding a system, spin-spin interactions could take on a number of different forms such as the Ising, flip-flop (or XX) and Heisenberg interactions in Eq. (6). For simplicity, we investigate the Ising interactionKrovi et al. (2007); Torrontegui and Kosloff (2016) in addition to the first form of interaction in Eq. (6).

The physical spin bath must contain a finite number of bath spins. In certain systems, such as electrically-gated quantum dotsHsieh et al. (2012) in GaAs, there could be as many as $10^{5}-10^{6}$ nuclear spins within the quantum-dot confining potential. While only a small fraction of bath spins are strongly coupled to the system; it is often possible to make semi-classical approximations to simplify the calculations. On the other hand, for NV centersGao et al. (2015) and related system, the relevant spin bath contains only $10^{1}-10^{2}$ spins. The bath could be potentially too small for semi-classical approximation and too large for a full dynamical simulations. Although we have seen impressive advances in simulation methodsWang and Shao (2012); Dobrovitski and De Raedt (2003); Al-Hassanieh et al. (2006); Stanek et al. (2013); Gelman et al. (2004) for large spin systems in the last decade, there still exists rooms for improvement.

II.3 Spin Bath Parameters

A finite-size bath model is fully characterized by pairs of parameters, $\{\omega_{k},g_{k}\}$. In modelling physical spin system, these parameters are often randomly drawn from narrow probability distributions as done and justified in earlier worksShenvi et al. (2005); Schliemann et al. (2003); Zhang et al. (2006). In particular, we will sample both $\omega_{k}$ and $g_{k}$ from the uniform distributions.

When addressing spin bath as a representation for anharmonic condensed phase environment, we consider an alternative assignment of parameters, $\{\omega_{k},g_{k}\}$, by discretizing an Ohmic spectral density, $J(\omega)\propto\omega\exp(-\omega/\omega_{c})$ according to the scheme given in Ref. Makri, 1999. In the thermodynamic limit, it has already been shownMakri (1999) that the spin bath can be exactly mapped onto a bosonic one with a temperature-dependent spectral density

[TABLE]

Our present focus is to investigate the nonlinear effects beyond the effective spectral density prescription.

III Methodology

III.1 Stochastic decoupling of many-body quantum dynamics

We now present an approach to systematically incorporate the non-linear bath effects into the HEOM framework through a stochastic calculus based derivation. Given the model described earlier, the exact quantum dynamics of the composite system (central spin and the bath) can be cast into a set of coupled stochastic differential equations,

[TABLE]

where $\tilde{\rho}_{s/B}(t)$ refers to the stochastically evolved density matrices in the presence of the white noises, implied by the Wiener differential increments $dW$ and $dV$, and the bath-induced stochastic field acting on the system,

[TABLE]

The reduced quantum dynamics of the central spin is recovered after averaging $\tilde{\rho}_{s}(t)$ over different noise realizations in the above equations,

[TABLE]

As implied in Eq. (13), the dissipative effects of the bath are completely captured by the interplay of the stochastic field $\mathcal{B}(t)$ and the white noises. To determined the stochastic field, one can use Eq. (III.1) to derive a closed form expression for $\mathcal{B}(t)$ in the case of bosonic bath,

[TABLE]

where $\alpha(t)=\int^{\infty}_{0}d\omega J(\omega)\left(\coth(\beta\omega/2)\cos(\omega t)-i\sin(\omega t)\right)$ is the two-time correlation functions. For non-Gaussian bath models, such as the spin bath, the stochastic field is determined by multi-time correlation functions as follows,

[TABLE]

where the definition on correlation functions $\Phi_{n,m}(t,t_{1},\dots,t_{n-1})$ are delegated to the App. A. The equation (18) assumes the odd-time correlation functions vanish with respect to the initial thermal equilibrium state. For now, we simply note that there are $2^{n-1}$ possible $n$-time correlation functions. The second subscript $m=1,\dots 2^{n-1}$ labels these functions. Not all of the $n$-time correlation functions are independent; every correlation function and its complex conjugate version are counted separately in this case. The motivation to distinguish the correlation functions is actually to differentiate all possible sequence of multiple integrations of noise variables associated with the $n$-time correlation functions as shown in Eq. (18). In this equation, we have explicitly expanded this formal expression up to the fourth-order correlation functions.

Substituting Eq. (18) into Eq. (13), one then derives a closed stochastic differential equation for the central spin. Direct stochastic simulation schemes so far have only been proposed for the Gaussian baths when the stochastic field is exactly defined by the second cumulant as in Eq. (17). In this case, the stochastic field can be combined with the white noises to define color noises with statistical properties specified by the bath’s two-time correlation function. When higher order cumulant terms are needed to properly characterize $\mathcal{B}(t)$, a direct stochastic simulation becomes significantly more complicated. In this study, we choose to convert Eqs. (13) and (18) to a hierarchy of deterministic equations involving auxiliary density matrices.

III.2 Generalized hierarchical equations of motion

To derive the hierarchical equation, we begin with Eq. (13) and take a formal ensemble average of the noises to get

[TABLE]

To arrive at the equation above, we invoke the relation of Eq. (16) and the fact $\mathcal{E}(dW)=\mathcal{E}(dV)=0$. This deterministic equation now involves an auxiliary density matrix $\mathcal{E}(\mathcal{B}(t)\tilde{\rho}_{s}(t))$ that needs to be solved too. Working out the equation of motion for the auxiliary density matrix (ADM), one is then required to define additional ADMs and a hierarchy forms.

Following a recently proposed scheme, we introduce a complete set of orthonormal functions $\{\phi_{j}(t)\}$ and express all the multi-time correlation functions as

[TABLE]

where $\boldsymbol{j}=(j_{1},\dots j_{n})$. Due to the completeness, one can also cast the derivatives of the basis functions in the form,

[TABLE]

Next we define the cumulant block matrices

[TABLE]

where each composed of $2^{n+1}$ row vectors with indefinite size and the [math] in each row implies all zeros beyond this point. For instance, $A_{1}$ has two row vectors while $A_{2}$ has four row vectors etc. The $m$-th row vector of matrix $A_{n}$ contains matrix elements denoted by $(a^{n}_{m\mathbf{j}_{1}},a^{n}_{m\mathbf{j}_{2}},\dots a^{n}_{m\mathbf{j}_{k}})$. Each of this matrix element can be further interpreted by

[TABLE]

where $dU(s_{j})$ can be either a $dW(s_{j})$ or $dV(s_{j})$ stochastic variable depending on index $m$. With these new notations, the multi-time correlation functions in Eq. (18) are now concisely encoded by

[TABLE]

Now we introduce a set of ADM’s

[TABLE]

which implies the noise average over a product of all non-zero elements of each matrix $\mathbf{A}_{i}$ with the stochastically evolved reduced density matrix of the central spin. The desired reduced density matrix would correspond to the ADM in which all matrices are empty. Furthermore, the very first ADM we discuss in Eq. (19) can be cast as

[TABLE]

where each ADM, $\rho^{\dots\left[\mathbf{A}_{n}\right]\dots}$, on the RHS of the equation carries only one non-trivial matrix element $a^{n}_{m,\boldsymbol{j}}$ in $\mathbf{A}_{n}$. Finally, The hierarchical equations of motion for all ADMs can now be put in the following form,

[TABLE]

In this equation, we introduce a few compact notations that we now explain. We use $\left[\mathbf{A}_{n}\pm(m,\boldsymbol{j})\right]$ to mean adding or removing an element $a^{n}_{m\boldsymbol{j}}$ to the $m$-th row. We also use $\left[a^{n}_{m\mathbf{j}}\rightarrow a^{n}_{m\mathbf{j^{\prime}}}\right]$ to denote a replacement of an element in the $m$-th row of $\mathbf{A}_{n}$. On the second line, we specify an element in a lower matrix given by $(m^{\prime},\mathbf{j}_{1})$. The variable $\mathbf{j}_{1}$ implies removing the first element of the $\mathbf{j}$ array and the associated index $m^{\prime}$ is determined by removing the first stochastic integral in Eq. (III.2). After the first term on the RHS of Eq. (III.2), we only explicitly show the matrices $\mathbf{A}_{n}$ affected in each term of the equation.

IV Results and Discussions

In the following numerical examples, we investigate the dynamics of a central spin coupled to a spin bath, Eq. (1). First two forms of system-bath interactions in Eq. (6) will be addressed. We use the gHEOM to simulate dynamics and adopt the Chebyshev polynomials as the functional basis to interpolate high-dimensional multi-time correlation functions. We note that there exists simpler choicesTang et al. (2016) of functional basis when only the two-time correlation functions are needed as in the bosonic bath models. We consider the standard initial conditions,

[TABLE]

where the bath density matrix is simply a tensor product of the thermal equilibrium state for each individual mode.

Through the examples in this section, we investigate whether it is generally a valid idea to map a spin bath model onto an effective bosonic one, such as obtained through a second order cumulant expansion of the spin bath’s influence functional. The gHEOM presented in Sec. III.2 will be used to quantify the contributions of higher order cumulant corrections to the quantum dynamics of the central spin.

IV.1 Pure dephasing

We first analyze a pure dephasing modelBhaktavatsala Rao and Kurizki (2011), i.e. $\Delta=0$ in Eq. (1) and adopt the first interaction $\hat{H}_{\text{int}}=\hat{\sigma}^{z}_{0}\sum_{k}g_{k}\hat{\sigma}_{k}^{x}$ in Eq. (6). This case can be analytically solved and provides insights into the higher order response functions of the bath. The coherence of the central spin can be expressed as

[TABLE]

where the decoherence factor $\Gamma(t)$ reads,

[TABLE]

with $\hat{H}^{(k)}_{\pm}=(\omega/2)\hat{\sigma}^{z}_{k}\pm g_{k}\hat{\sigma}^{x}_{k}$ and $\Omega_{k}=\omega_{k}\sqrt{1+(4g_{k}^{2}/\omega_{k}^{2})}$. On the last line, we expand $\Gamma(t)$ to get the two leading contributions with respect to $\lambda_{k}=4g^{2}_{k}/\omega_{k}^{2}$. These two terms correspond to the second order and fourth order cumulant expansion of the influence functional for this particular model, respectively.

Even with this simple case, one can draw important remarks regarding spin bath mediated decoherence. First of all, the exact result in Eq. (33) implies the perturbations coming from a specific spin bath mode are modulated with an interaction renormalized frequency $\Omega_{k}$ as opposed to the bare frequency $\omega_{k}$, which is the way bosonic modes perturb a system through a linear coupling. The origin of this interaction dressed $\Omega_{k}$ can be understood by inspecting the time evolution of Pauli matrices associated with individual spin modes,

[TABLE]

By converting the raising and lowering Pauli matrices into the $x$ and $y$ Pauli matrices, it is obvious that different components of the bath spin get coupled together via the system-bath interaction in a non-trivial way and renormalize the frequency at which a spin precesses. There is no such coupling of the internal structure for harmonic oscillators due to the fundamental differences in the commutation properties of bosonic creation/annihilation operators and Pauli matrices for spins.

Secondly, when $\lambda_{k}\gg 1$ then $\Omega_{k}$ is significantly shifted from the bare frequency $\omega_{k}$. The density of states of the bath will be dramatically re-organized with respect to $\Omega_{k}$. Any methods expanding around $\omega_{k}$ will be difficult to provide accruate results when system-bath coupling is strong and/or when $\omega_{k}$ is small. A striking example would be a single-frequency bath, in which all modes possess the identical energy scale $\omega_{k}=\omega$ and coupled non-uniformly to the system. For a bosonic bath as well as the second-order cumulant expansion for a spin bath, there is essentially no dephasing. According to the second-order result in Eq. (33), the coherence of the central spin is periodically recovered at time points $\omega t=n2\pi$ with $n$ an integer. However, for an exact treatment of the spin bath, non-trivial dephasing happens as long as not all coupling coefficients $g_{k}$ are identical. Due to the system-bath interaction, the single-frequency distribution of $\omega$ can be broadened to a finite bandwidth corresponding to the dressed $\Omega_{k}$. This analysis is confirmed in Fig. 1, where the exact result (green) undergo an irrversible decay but not for the second-order result (red). In the figure, we also consider a modified second order cumulant result (blue dotted) which expands the decoherence function $\Gamma(t)$ with respect to $\Omega_{k}$ instead of $\omega_{k}$. As shown, this modified expansion works extremely well. In general, this dressing of $\Omega_{k}$ implies a faster dephasing rate is expected from a spin bath when compared to a similar bosonic bath, i.e. a bath of harmonic oscillators sharing the same set of $\{\omega_{k},g_{k}\}$.

Finally, the temperature independence of $\Gamma(t)$ in Eq. (33) is another distinguishing property to set apart spin bath from what has been known for the bosonic bath models. This temperature-independent dephasing can already be inferred from the expression of the effective spectral density in Eq. (12), and is further confirmed to hold beyond the linear response regime in this pure dephasing model as the temperature factors is missing in the exact expression in Eq. (33).

To estimate the contributions of higher order cumulant corrections (HOCCs) to the central spin deocherence, we note the fourth order cumulants in the last line of Eq. (33) can become dominant under two conditions: (1) the perturbation parameters satisfy $\lambda_{k}>1$ and/or (2) when $t>1/(\lambda_{k}^{2}\omega_{k})$ such that the terms linearly proportional to $\omega_{k}t$ dominates the second order cumulant. The second condition implies a potential linear time $t$ divergence. This instability is an artifact of cumulant expansion and can be removed by introducing higher order cumulants. When the bath parameters $\{\omega_{k},g_{k}\}$ are obtained by discretizing a continuous spectral density as discussed in Sec. II.3, all $\lambda_{k}\sim 1/N_{B}$ can be made arbitrarily small when sufficiently large number of spin modes are used. Hence, the artificial divergence due to unstable part of the fourth-order cumulants can often be suppressed within the typical time domain for simulating condensed-phase dynamics. However, even if just a few modes, satisfying $\lambda_{k}\geq 1$, could potentially contribute immensely to the overall HOCCs because of the logarithmic form for each spin’s contribution to the exact expression for $\Gamma(t)$ in Eq. (33).

As mentioned in the case of physical spin modelsProkofiev and Stamp, P. C. E. (2000), there is no reason that $\lambda_{k}$ should be related to the number of bath spins. Hence, the effects of HOCCs can become noticeable if not all $\lambda_{k}$ are sufficiently small within the simulation time window to suppress the divergence associated with the unstable part of higher order cumulants. In Fig. 2, we look at the dephasing rate, $\Gamma(t)/t$ for two different cases to analyze HOCCs. The initial condition of the central spin is taken to be the pure state $\left|\psi\right\rangle=\frac{1}{\sqrt{2}}\left(\left|\uparrow\right\rangle+\left|\downarrow\right\rangle\right)$. This calculation also serves as a benchmark to validate that the gHEOM correctly resolves the second and fourth cumulant contributions when compared to the exact expansions in Eq. (33). The parameters $\epsilon=2$ and $\Delta=0$ are used for the system Hamiltonian. We assign random samples of $\{\omega_{k},g_{k}\}$ from two uniform distributions centered on $\omega_{o}$ and $g_{o}$, respectively, with details given in the figure caption. In panel (a), the fourth order results can sufficiently reproduce the exact rate and it starts to deviate from the linear response rate around $t\sim 3$. In panel (b), by further increasing the average coupling $g_{0}$, even the fourth order corrections starts to fall short of reproducing the exact rate around $t\sim 3$. For both cases considered in Fig. 2, $\lambda_{k}<1$ hold for all bath modes, which ensure all perturbative expansion parameters $\lambda_{k}$ are well-behaved in the short-time limit. It is clear that the HOCCs becomes important to simulate the system relaxation when the spin bath parameters do not satisfy the scaling relation $\lambda_{k}\sim 1/N_{B}$.

IV.2 Anharmonic Condensed-Phase evironment

We consider the same model as in Sec. IV.1 but with $\Delta\neq 0$ in Eq. (1). The central spin could suffer relaxation due to interaction with the bath in this case. The parameters,$\{\omega_{k},g_{k}\}$, are assigned by discretizing an Ohmic spectral density. This finite-size restriction is a necessity to observe any deviations from linear response results as explained in the introduction. Furthermore, the discretization allows us to numerically compute the multi-time correlation functions, needed for the gHEOM calculations, by summing over contributions from each bath spins.

One can estimate the leading order corrections beyond the linear response approximation by a perturbative expansion with respect to $\lambda_{k}$ as done in the previous section. We first analyze the population dynamics in the Markovian limit. We adapt the NIBA (Non-Interacting Blip Approximation) equation to the spin bath modelSegal (2014) with a symmetric system Hamiltonian, $H_{s}=\frac{\Delta}{2}\sigma^{x}_{0}$. We find

[TABLE]

where the $Q_{R}$ and $Q_{I}$ functions read

[TABLE]

with $\lambda_{k}\equiv(4g_{k}/\omega_{k})^{2}$. The equation (35) is a second-order expansion (with respect to the off-diagonal element, $\Delta$) of the memory kernel. The functions, $Q_{R}(t)$ and $Q_{I}(t)$, in Eq. (IV.2) are further expanded up to $\lambda_{k}^{2}$. If one only retains the first term, proportional to $\lambda_{k}$, then $Q_{R}(t)$ and $Q_{I}(t)$ reduce to the standard NIBA expressions for an effective bosonic bath with a temperature dependent spectral density, Eq. (12).

To finish the Markovian approximation, we replace $<\sigma^{z}_{0}(s)>$ with $<\sigma^{z}_{0}(t)>$ on the RHS of Eq. (35), extend the integration limit to infinity in both directions, and perform a short-time expansions of $Q_{R}$ and $Q_{I}$ to keep terms up to $\sim t^{2}$. One then integrates out the memory kernel in Eq. (35) and obtains a simple rate equation with the Fermi golden rule rate given by

[TABLE]

where $\xi=a_{2}/a_{1}$, $\zeta=b_{2}/b_{1}$, $a_{1}=\sum_{k}\lambda_{k}\omega_{k}^{2}/2$, $a_{2}=\sum_{k}\lambda_{k}^{2}\tanh(\beta\omega_{k}/2)\omega_{k}^{2}/2$, $b_{1}=\sum_{k}\lambda_{k}\omega_{k}\tanh(\beta\omega_{k}/2)$ and $b_{2}=\sum_{k}\lambda_{k}^{2}\omega_{k}\tanh(\beta\omega_{k}/2)/2$. As $\lambda_{k}\rightarrow 0$, the rate approaches to the linear response / effective bosonic bath results: $k\rightarrow k_{\text{lin}}$. The expression on the second and third line constitutes a good approximation when both $\xi$ and $\zeta$ are small. In particular, the last exponential factor isolates the leading-order correction to the rate constant, $\eta=\exp\left(-\frac{1}{2}\xi\right)$. Whenever the scaling $\lambda_{k}\sim 1/N_{B}$ is imposed, $\xi$ will scale as $1/N_{B}$ too and $\eta$ will be exponenially suppressed. Furthermore, we note that the leading-order correction reduces the relaxation rate at low temperature and gradually converge to the linear response result $k_{\text{lin}}$ with increasing temperature due to the factor $\tanh(\beta\omega_{k}/2)$ contained in $a_{2}$ variable.

Next, we investigate numerically the convergence of a spin bath to the linear response results within the gHEOM calculations. We use the following parameters $\Delta=1$ and $\epsilon=0$ for the system Hamiltonian and the spectral density parameters: $\omega_{c}=\Delta$, $\beta\Delta=2$, $\alpha=2.3$. For Ohmic spectral density, coupling coefficients satisfy $g_{k}^{2}\propto\omega_{k}$ and the perturbative parameter $\lambda_{k}\propto 1/\omega_{k}$. Therefore, the low frequency bath modes are more severely influenced by the system-bath interactions with $\Omega_{k}=\omega_{k}\sqrt{1+\lambda_{k}}$ shifted further from $\omega_{k}$. To investigate deviations from the linear response results, the highest frequency we consider is $\omega_{\text{max}}=2\omega_{c}$ in the discretization. The convergence of dynamics is demonstrated in Figs. 3a and 3b, the population and the coherences of the central spin are plotted up to $\Delta t=3.5$, respectively, for a total number of $35$, $70$, $105$, $210$ and $500$ spins including up to the fourth-order cumulant expansions. As the number of bath spins increase, the results converge smoothly to those of the condensed phase environment in the thermodynamic limit. At $N_{B}=500$, the fourth order results already converge with the second order results.

Next, we investigate the temperature dependence of HOCCs. The same set of Hamiltonain and bath’s spectral density parameters as above is used except the temperature will be varied for analyses. We use a set of $N_{B}=35$ spins for illustrations. In many earlier studies, distinctive properties of a spin bath (in comparison to a bosonic counterpart) are found to be temperature-related and are attributed to the temperature-dependent spectral density, Eq. (12). However, restricting the discussions to the effective spectral density imply the comparisons focused on the linear response limit. We would like to further analyze the contributions of HOCCs to these temperature-dependent effects. In Fig. 4, we compare the second-order (dashed curves) and fourth-order (solid curves) results to inspect the contribution of the HOCCs. In short, the HOCCs become more pronounced when the temperature is lowered and the divergence between second-order and four-order results increase. The numerical results is also consistent with the earlier conclusion drawn from the rate expression, Eq. (37). Similar observationsMakri (1999); Wang and Shao (2012) have been reported in the literature where it was found that more number of discretized bath spins are needed to reach the linear response limit at low temperatures. Despite the simple argument that the linear response of a spin and a harmonic oscillator converge in the zero temperature limit, the two bath models actually do not converge except in the cases of large bath size. This is because the higher order response functions for spins become more prominent in the low temperature regime. The primary reason is due to the way temperature enters the correlation functions as $\tanh(\beta\omega_{k}/2)$.

IV.3 Ising Spin Bath

Finally, we consider another system-bath interaction, $\hat{H}_{\text{int}}=\hat{\sigma}_{0}^{z}\sum_{k}g_{k}\hat{\sigma}_{k}^{z}$. The Ising spin-spin interaction prohibits bath spins to be flipped but still entangles system and bath. While this spin bath model is appropriate in certain quantum computing contextsKrovi et al. (2007); Torrontegui and Kosloff (2016), it does not relate to a condensed phase environment. Hence, we will follow the physical spin bath approach to sample $\omega_{k}$ and $g_{k}$ from uniform distributions.

We first consider a pure dephasing case with $\Delta=0$. The coherence of the central spin can be cast in the general form of Eq. (32) but with a different decoherence function,

[TABLE]

where $\gamma_{k}=\tanh(\beta\omega_{k}/2)$. Unlike the previouse pure-dephasing model in Sec. IV.1, $\Gamma(t)$ acquires a temperature dependence through $\gamma_{k}$. Since $\hat{H}_{B}$ and $\hat{H}_{\text{int}}$ commutes, the bath Hamiltonian can be removed from the dynamical equation in a rotated frame and no dressed bath frequencies $\Omega_{k}$ appear as in Sec. IV.1. For the Ising spin bath, the bath frequencies $\omega_{k}$ enter the decoherence function through $\gamma_{k}$, reflecting the thermal equilibrium initial condition $\rho^{eq}_{B}$.

According to Eq. (38), $\gamma_{k}$ modulate the magnitude of $\Gamma(t)$ and the Ising spin bath becomes less efficient at interacting with the system in the high temperature limit and / or the low-frequency limit when $\gamma_{k}\ll 1$. The system-bath coupling coefficients $g_{k}$ determine the oscillatory behaviors of $\Gamma(t)$ in time domain. When $g_{k}$ are narrowly distributed around a mean value $g_{o}$, one expects a partially periodic recurrence of $\Gamma(t)$.

Now we investigate effects of HOCCs. Two points make the Ising spin model an interesting case to analyze. First, an expansion of Eq. (38) with respect to $g_{k}$ reveals that the even cumulants do not contribute to the imaginary component of $\Gamma(t)$, which drives a rotation of the central spin on the Bloch sphere. A second order cumulant expansion will completely miss this rotation. This point is illustrated in Fig. 5 where the central spin’s initial condition is taken to be $\left|\psi_{s}(0)\right\rangle=(\left|\uparrow\right\rangle+\left|\downarrow\right\rangle)/\sqrt{2}$. The second-order result (red curve) reproduces the real part of the coherence, $\langle\uparrow|\rho_{s}(t)|\downarrow\rangle$, but completely misses the growth of the imaginary component. Once we incorporate the third and fourth cumulants in the calculation, the result (green curve) converges better to the exact one in Fig. 5b. Secondly, take $B=\sum_{k}g_{k}\hat{\sigma}_{k}^{z}$ as the bath part of $\hat{H}_{\text{int}}$ and it is easy to verify the multi-time correlation functions, such as $\text{Tr}\left\{\rho^{eq}_{B}B(t_{1})B(t_{2})\dots B(t_{n})\right\}$, are time invariant (i.e. the memory kernel of the bath do not decay in time). The highly non-Markovian nature of the Ising spin bath can give rise to non-trivial steady states in the long time limit. In Fig. 6, we consider another case and find the second-order result suggests a fully decayed coherence while the exact result shows a persistent oscillation. Although the higher-order result (green curve) can better capture the on-going oscillating behavior in the transient regime, the artificial divergence of cumulant expansion (explained in Sec. IV.1) suggests more cumulant terms should be included within the simulated time domain.

Going beyond the pure dephasing case, we restore the off-diagonal coupling $\Delta$. We examine both the coherence and population dynamics in Fig. 7. Similar to pure dephasing case, one expects a slow decoherence when $g_{k}$ are narrowly distributed. For both coherence and population dynamics, we identify the clear insufficiency of second-order results and higher-order cumulants again helps to restore the coherence and the beatings of the population dynamics in the transient regime.

V Conclusions

In summary, we use our recently proposed gHEOM to investigate the effects of HOCCs on the quantum dissipations induced by finite-size spin bath models. The gHEOM can systematically incorporate the higher order cumulants of the bath’s influence functional into calculations. The controlled access to non-Gaussian effects of the bath allows us to assess the sufficiency of a linear response approximation. Besides the spin baths, the methodology can be similarly applied towards other types of anharmonic environments. However, due to the prohibitive numerical resources required to accurately characterize the high-dimensional multi-time correlations functions, the best usage of this method is to combine it with transfer tensor method (TTM) proposed by one of us. One can use the gHEOM to quantitatively capture the exact short-time dissipative dynamics embedded in a non-Gaussian bath. These short-time results are then fed to the TTM method to reproduce the correct memory kernel of the environment and allow an efficient and stable long-time simulation of dissipative dynamics.

Through the analyses done in Sec. IV.2, we find the linear response approximation provides a highly efficient and accurate result for a finite spin bath over a wide range of parameters. This is mainly because the next leading order correction scale as $1/N_{B}$ in the cumulant expansion. We present one “extreme” result for a relatively slow Ohmic bath in order to observe appreciable corrections coming from the higher order cumulant terms in the short time limit. Even in this case, the higher order effects still vanish when $N_{B}=500$. Although, the low-temperature condition should exacerbate the discrepancy between exact and linear-response results, we find the actual effects rather minimal in the short-time limit. Considering the significant numerical costs to access higher order cumulants, it certainly make linear response approximation a highly appealing option in dealing with a spin-based condensed phase environment. In App. B, we further investigate the differences between a spin and bosonic bath in the linear response limit. We confirm the lack of appreciable temperature dependence on the dissipative dynamics and the emergence of negative differential thermal conductance are two robust physical signatures to distinguish a spin bath from a corresponding bosonic one as explained in the appendix.

It is much simpler to devise numerical examples in which the higher order cumulants play critical roles in a physical spin bath model. The most critical factor is the probability distributions for $\{\omega_{k},g_{k}\}$. Narrow distributions will make the spin bath model more difficult for linear-response approximations, and it is likely to have such narrow distributions in real spin-based environments. In such cases, linear-response results could deviate extremely from the exact results such as shown in Fig. 1 and Fig. 6. Secondly, for highly symmetric spin-spin interaction such as the Ising Hamiltonian considered in Sec. IV.3, the second order cumulants fail to generate a rotation of the spin state, which could only be accounted by the odd-order cumulants. Finally, the physical spin bath could be difficult to handle due to the possibility of extreme non-Markovianity. In the Ising Hamiltonian example, we see the extreme case of having all bath’s multi-time correlation functions to be time-invariant and we need to expand deep down the hierarchy to obtain converged results. The highly non-Markovian nature of spin bath is not a rare exception. In addition to Ising Hamiltonian, the flip-flop and Heisenberg Hamiltonian in combination with narrow distribution of $\{\omega_{k},g_{k}\}$ can also result in highly symmetric systems (central spin plus the bath) with a rich set of non-Markovian and persistent dynamics. With the gHEOM method, we can systematically incorporate higher order cumulants to improve the simulation results for physical spin bath models.

Acknowledgements.

C.H. acknolwedges support from the SUTD-MIT program. J.C. is supported by NSF (grant no. CHE-1112825) and SMART.

Appendix A Stochastic mean field and multi-time correlation functions

From Eq. (15), it is clear that equation of motions for $\mathcal{B}(t)$ can be obtained from the stochastic dynamical equations for the bath density matrix in Eq. (III.1). More specifically, the time evolution of the stochastic field $\mathcal{B}=\sum_{k}g_{k}(\langle b^{{\dagger}}_{k}\rangle+\langle b_{k}\rangle)$ is jointly determined by

[TABLE]

The expectation values in Eqs. (39)-(40) are taken with respect to the stochastically evolved $\tilde{\rho}_{k}(t)$. The generalized cumulants above are defined as

[TABLE]

In this case of bosonic bath models, it is straightforward to show that the time derivatives of $\mathcal{G}^{\alpha_{1}\alpha_{2}}_{k}$ vanish exactly. Hence, the second order cumulants are determined by the thermal equilibrium conditions of the initial states. Immediately, one can identify the relevant quantity $\mathcal{G}^{+-}_{k}=n_{B}(\omega_{k})$, the Bose-Einstein distribution for the thermal state of the bath. The time invariance of the second order cumulants make Eq. (39)-(40) amenable to deriving a closed form solution. On the other hand, for non-Gaussian bath such as the spin bath, the second cumulants are not time invariant. One way to determine their time evolution is to work out their equations of motion by iteratively applying Eq. (III.1). It is straightforward to show these equations couples different orders of generalized cumulants,

[TABLE]

where $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2}\dots\alpha_{n})$ specifies a sequence of raising and lowering spin operators that constitute this particular $n$-th order cumulant and $|\boldsymbol{\alpha}|=\sum_{i}\alpha_{i}$ with $\alpha_{i}=\pm 1$ depending on whether it refers to a raising (+) or lowering (-) operator, respectively. We use $[\boldsymbol{\alpha},\pm]\equiv(\alpha_{1},\dots\alpha_{n},\pm)$ to denote an $n+1$-th cumulant obtained by appending a spin operator to $\boldsymbol{\alpha}$. A similar definition is implied for $[\pm,\boldsymbol{\alpha}]$. More specifically, these cumulants are defined via an inductive relation that we explicitly demonstrate with an example to obtain a third-order cumulant starting from a second-order one given in Eq. (41),

[TABLE]

The key step in this inductive procedure is to insert an operator identity $b^{\pm}-\langle b^{\pm}\rangle$ at the end of each expectation bracket defining the $n$-th cumulant. If a term is composed of m expectation brackets, then this insertion should apply to one bracket at a time and generate m terms for the $n+1$-th cumulant. Similarly, we get $\mathcal{G}^{[\pm,\boldsymbol{\alpha}]}$ by inserting the same operator identity to the beginning of each expectation bracket of $\mathcal{G}^{\boldsymbol{\alpha}}_{k}$.

For the spin bath, these higher order cumulants persists up to all orders. In any calculations, one should certainly truncate the cumulants at a specific $k$-th order by imposing the time invariance, $G^{\boldsymbol{\alpha}}_{k}(t)=G^{\boldsymbol{\alpha}}_{k}(0)$ and evaluate the lower order cumulants by recursively integrating Eq. (42). Through this simple prescription, one derives Eq. (18).

Appendix B Physical signatures of the spin bath models in the thermodynamic limit

In the main text, we focus predominantly on the higher order corrections to the quantum dissipations induced by a spin bath. Now we turn attention to the linear response limit, and we look for physical signatures that can distinguish the spin and bosonic bath models . The condensed-phase spin bath (with practically infinite number of modes) can be rigorously mapped onto an effective bosonic bath with a temperature-dependent spectral density, Eq. (12). To facilitate the calculations in this appendix, we should take the spin bath as an effective bosonic model and adopt the superohmic spectral density for convenience. The results below compliments those of earlier studiesSegal (2014); Wang and Shao (2012); Gelman et al. (2004); Zhang et al. (2015); Lü and Zheng (2009) on the same subject.

B.1 Electronic Coherence of a Two-Level System

We first compare the dynamics of a dimer coupled to (1) a bosonic bath and (2) a spin bath. The spectral densities for the two models read

[TABLE]

where the subscript $b/s$ denotes the bosonic bath and the spin (effective bosonic) bath model, respectively. In the study of FMO-like molecular systems, surprisingly long coherence times were discovered and successfully explained through an enhanced NIBA formalismU (2012) which takes into account the first order blip interactions. By considering the same dimer system and a similar set of experimentally relevant parameters as in Ref. Pachón and Brumer, 2011, we investigate how much the coherent population dynamics of an effective dimer will change when the environment is replace by a spin bath with an identical set of parameters as the original harmonics-oscillator based condensed phase. The enhanced NIBA formalism is also adopted here as it provides sufficiently accurate results over a long time span in the parameter regime considered here. The parameters are $K=0.16$, $\epsilon/\Delta=0.6$, $\omega_{c}/\Delta=2.0$ and $\Delta\approx 106.2\text{cm}^{-1}$.

From Fig. 8, the dimer’s population dynamics in the presence of the bosonic (panel a) and the spin (panel b) bath at two different temperatures: $T=290K$ and $75K$ are presented. The significant temperature-dependent relaxation is observed in the standard bosonic bath; while the almost temperature independent relaxation is found in the spin bath model. These results confirm that the atypical temperature dependence of a spin bath induced relaxation can be quite pronounced and easily detected in the experimentally accessible regime.

B.2 Energy transport through a non-equilibrium junction

We next explore additional features of a spin-based environment in the non-equilibrium situations. We study the energy transport across a molecular junction connected to two heat baths composed of non-interacting spins. The extended double-bath model should read,

[TABLE]

where $\hat{H}_{s}$ and $\hat{H}_{B,\mu}$ are the standard system and bath Hamiltonian with $\mu=L,R$ to denote the two bath at left or right end. Recent studySegal (2014) tried to model anharmonic junctions with the spin bath and hinted several qualitative differences in the transport phenomena. In this work, we adopt a non-equilibrium Polaron-transformed Redfield equation (NE-PTRE) in conjunction with the full counting statistics to compute the steady-state energy transfer through the junction. In Ref.Wang et al., 2015, it was demonstrated that NE-PTRE can be reliably used to calculate energy currents (through a molecular junction coupled to bosonic baths) from the weak to the strong system-bath coupling regimes as the corresponding analytical expressions for the energy current reduced elegantly to either the standard Redfield (weak coupling) or NIBA (strong coupling) results in the appropriate limits. In this study we apply the NE-PTRE formalism to calculate and compare the energy current through the junction while contacted by (a) two spin baths and (b) two bosonic baths held at different temperatures. Similar to App. B.1, the actual calculations below will treat the spin bath as an effective bosonic bath. The same superohmic spectral density, Eq. (B.1), is adopted here. The parameters are $\Delta=1$, $\epsilon=10$, $\omega_{c}=10$ and $K=3.5$. The left and right bath share an identical set of parameters except the temperature, and the fast bath (or scaling) limit with $\omega_{c}\gg 1$ is imposed in both bath models.

In Fig. 9, the heat transfer exhibits a negative differential thermal conductance (NDTC) for spin baths but not for bosonic baths. The temperature, $k_{B}T_{\alpha}$ are measured in terms of $\Delta$ with $k_{B}=1$ in the present case. We remark that the NDTC (for bosonic bath models) reported in some earlier studies are found to be an artifact of the Marcus approximation (over a wide range of parameter regimes) as clarified by the NE-PTRE method reported in Ref. Wang et al., 2015. In Fig. 9, the NE-PTRE predicts the NDTC for the spin bath model but not for the bosonic bath model. In a separate work (not shown here), we also investigate a single-frequency spin bath and compute the energy current by using the NIBA method and find the same NDTC appearing in the strong coupling regime. This helps to validate that the present result is not specific to a particular spectral density. This qualitative difference between the bosonic and spin baths results make the NDTC another strong physical signature to distinguish the two bath models.