A Unified Stochastic Formulation of Dissipative Quantum Dynamics. I. Generalized Hierarchical Equations

TL;DR

This paper develops a unified stochastic framework to exactly simulate dissipative quantum dynamics across bosonic, fermionic, and spin environments, extending hierarchical equations of motion to non-Gaussian baths.

Contribution

It introduces a generalized stochastic Liouville equation and hierarchical equations that unify and extend existing methods for simulating open quantum systems with various bath types.

Findings

Exact simulation of quantum dissipation for different bath models.

Reduction to existing hierarchical equations for bosonic baths.

Efficient simulation methods for fermionic and spin baths.

Abstract

We extend a standard stochastic theory to study open quantum systems coupled to generic quantum environments including the three fundamental classes of noninteracting particles: bosons, fermions and spins. In this unified stochastic approach, the generalized stochastic Liouville equation (SLE) formally captures the exact quantum dissipations when noise variables with appropriate statistics for different bath models are applied. Anharmonic effects of a non-Gaussian bath are precisely encoded in the bath multi-time correlation functions that noise variables have to satisfy. Staring from the SLE, we devise a family of generalized hierarchical equations by averaging out the noise variables and expand bath multi-time correlation functions in a complete basis of orthonormal functions. The general hiearchical equations constitute systems of linear equations that provide numerically exact simulations of quantum dynamics. For bosonic bath models, our general hierarchical equation of motion reduces exactly to an extended version of hierarchical equation of motion which allows efficient simulation for arbitrary spectral densities and temperature regimes. Similar efficiency and exibility can be achieved for the fermionic bath models within our formalism. The spin bath models can be simulated with two complementary approaches in the presetn formalism. (I) They can be viewed as an example of non-Gaussian bath models and be directly handled with the general hierarchical equation approach given their multi-time correlation functions. (II) Alterantively, each bath spin can be first mapped onto a pair of fermions and be treated as fermionic environments within the present formalism.

Figures

Equations

$$\hat{H}$$

$${\hat{H}}_s=\frac{{\omega }_0}{2}{\hat{\sigma }}_0^z+\frac{\Delta }{2}{\hat{\sigma }}_0^x,$$

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

$$\frac{\partial \hat{\rho }}{\partial t}=-i\left[\hat{H},\hat{\rho }\right].$$

$$d{\tilde{\rho }}_s=$$

$$-idt\left[{\hat{H}}_s,{\tilde{\rho }}_s\right]-\frac{i}{\sqrt{2}}\sum _{\alpha =\pm }{A}^{\alpha }{\tilde{\rho }}_{s}d{W}_{\alpha }+\frac{i}{\sqrt{2}}\sum _{\alpha =\pm }{\tilde{\rho }}_s{A}^{\alpha }d{V}_{\alpha },$$

$$d{\tilde{\rho }}_B=$$

$$-idt\left[{\hat{H}}_B,{\tilde{\rho }}_B\right]+\frac{1}{\sqrt{2}}\sum _{\alpha =\pm }d{W}_{\alpha }^{\ast }{B}^{\bar{\alpha }}{\tilde{\rho }}_B+\frac{1}{\sqrt{2}}\sum _{\alpha =\pm }d{V}_{\alpha }^{\ast }{\tilde{\rho }}_B{B}^{\bar{\alpha }},$$

$$d{\tilde{\rho }}_k=-idt\left[H_{b,k},{\tilde{\rho }}_k\right]$$

$$+\frac{1}{\sqrt{2}}\sum _{\alpha =\pm }{g}_{k}d{W}_{\alpha }^{\ast }{b}_k^{\bar{\alpha }}{\tilde{\rho }}_k+\frac{1}{\sqrt{2}}\sum _{\alpha =\pm }{g}_{k}d{V}_{\alpha }^{\ast }{\tilde{\rho }}_k{b}_k^{\bar{\alpha }},$$

$$\overline{{\mu }_{\alpha }(t)}=\overline{{\nu }_{\alpha }(t)}=0,$$

$$\overline{{\mu }_{\alpha }(t){\mu }_{{\alpha }^{\prime }}^{\ast }(t^{\prime })}=\overline{{\nu }_{\alpha }(t){\nu }_{{\alpha }^{\prime }}^{\ast }(t^{\prime })}=2{\delta }_{\alpha ,{\alpha }^{\prime }}\delta (t-t^{\prime }),$$

$$d\tilde{\rho }(t)=[d{\tilde{\rho }}_s(t)]{\tilde{\rho }}_B(t)+{\tilde{\rho }}_s(t)[d{\tilde{\rho }}_B(t)]+d{\tilde{\rho }}_s(t)d{\tilde{\rho }}_B(t),$$

$${\hat{\rho }}_s(t)=\overline{{\tilde{\rho }}_s(t){\text{Tr}}_B{\tilde{\rho }}_B(t)}.$$

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

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

$$d{\tilde{\rho }}_B=-idt\left[{\hat{H}}_B,{\tilde{\rho }}_B\right]+\frac{1}{\sqrt{2}}\sum _{\alpha }d{W}_{\bar{\alpha }}^{\ast }\left(B^{\alpha }\mp {\mathcal{B}}^{\alpha }(t)\right){\tilde{\rho }}_B$$

$$+\frac{1}{\sqrt{2}}\sum _{\alpha }d{V}_{\bar{\alpha }}^{\ast }{\tilde{\rho }}_B\left(B^{\alpha }-{\mathcal{B}}^{\alpha }(t)\right).$$

$$\displaystyle\mathcal{B}^{\alpha}(t)=\left\{\begin{array}[]{ll}\sum_{k}g_{k}\text{Tr}_{B}\left\{\tilde{\rho}_{B}(t)b^{\alpha}_{k}\right\},&\text{complex-valued}\\ \alpha\sum_{k}g_{k}\text{Tr}_{B}\left\{\tilde{\rho}_{B}(t)b^{\alpha}_{k}\right\},&\text{Grassmann-valued},\end{array}\right.$$

$$d⟨b_k^{†}⟩$$

$$d⟨b_k⟩$$

$${\mathcal{G}}_k^{{\alpha }_1{\alpha }_2}=⟨b_k^{{\alpha }_1}{b}_k^{{\alpha }_2}⟩\mp ⟨b_k^{{\alpha }_1}⟩⟨b_k^{{\alpha }_2}⟩.$$

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

$$\displaystyle\alpha(t)=\left\{\begin{array}[]{ll}\sum_{k}|g_{k}|^{2}\left(\cos(\omega_{k}t)\coth(\beta\omega_{k}/2)-i\sin(\omega_{k}t)\right),&\text{bosonic bath}\\ \sum_{k}|g_{k}|^{2}\left(\cos(\omega_{k}t)-i\sin(\omega_{k}t)\tanh(\beta\omega_{k}/2)\right),&\text{fermionic bath}\\ \end{array}\right.$$

$$d{\mathcal{G}}_k^{\mathbf{\alpha }}$$

$${\mathcal{G}}_k^{[({\alpha }_1,{\alpha }_2),\pm ]}=⟨b^{{\alpha }_1}{b}^{{\alpha }_2}(b^{\pm }-⟨b^{\pm }⟩)⟩+⟨b^{{\alpha }_1}(b^{\pm }-⟨b^{\pm }⟩)⟩⟨b^{{\alpha }_2}⟩+⟨b^{{\alpha }_1}⟩⟨b^{{\alpha }_2}(b^{\pm }-⟨b^{\pm }⟩)⟩.$$

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

$$\frac{d{\tilde{\rho }}_s}{dt}=-i\left[{\hat{H}}_s,{\tilde{\rho }}_s\right]\mp iA{\tilde{\rho }}_s(t)\xi (t)+i{\tilde{\rho }}_s(t)A\eta (t),$$

$$\xi (t)$$

$$\eta (t)$$

$$\overline{\xi (t)\xi (t^{\prime })}$$

$$\overline{\xi (t)\eta (t^{\prime })}$$

$$\overline{\eta (t)\eta (t^{\prime })}$$

$$\xi (t)$$

$$\eta (t)$$

$$\alpha (t-t^{\prime })$$

$${\alpha }_R(|t-t^{\prime }|)$$

Full text

I Introduction

Understanding 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 have been achieved within the context of a few prototypical modelsLeggett et al. (1987); Breuer and Petruccione (2002); Weiss (2012) such as Caldeira-Leggett model and spin-boson model. In most cases the environment is modeled as a bosonic bath, a set of non-interacting harmonic oscillators whose influences on the system is concisely encoded in a spectral density. The prevalent adoption of bosonic bath models is based on the arguments that knowing the linear response of an environment near equilibrium should be sufficient to predict the dissipative quantum dynamics of the system.

Despite many important advancements in the quantum dissipation theory have been made with the standard bosonic bath models in the past decades, more and more physical and chemical studies have suggested the essential roles that other bath models assume. We briefly summarize three scenarios below.

A standard bosonic bath model fails to predict the correct electron transfer rate in donor-acceptor complex strongly coupled to some low-frequency intramolecular modes. Some past attempts to model such an anharmonic, condensed phase environment include (a) using a bath of non-interacting MorseWang and Thoss (2007); López-López et al. (2011); Kryvohuz and Cao (2005); Wu and Cao (2001) or quartic oscillatorsand (b) mapping anharmonic environment onto effective harmonic modesCao and Voth (1995); Makri (1999) with a temperature-dependent spectral density. 2. 2.

Another prominent example is the fermonic bath model. Electronic transports through nanostructures, such as quantum dots or molecular junctions, involves particle exchange occurs across the system-bath boundary. Recent developments of several many-body physics and chemistry methods, such as the dynamical mean-field theoryKotliar et al. (2006) and the density matrix embedding theoryKnizia and Chan (2012), reformulate the original problem in such a way that a crucial part of the methods is to solve an open quantum impurity model embedded in a fermionic environment. 3. 3.

The spin (two-level system) bath models have also received increased attention over the years due to ongoing interests in developing various solid-state quantum technologiesHsieh et al. (2012) under the ultralow temperature when the phonon or vibrational modes are frozen and coupling to other physical spins (such as nuclear spins carried by the lattice atoms), impurities or defects in the host material emerge as the dominant channels of decoherence.

Both bosonic and fermionic environments are Gaussian baths, which can be exactly treated by the linear responseMakri (1999) in the path integral formalism. For the non-Gaussian baths, attaining numerically exact open quantum dynamics would require either access to higher order response function of the bath in terms of its multi-time correlation functions or explicit dynamical treatments of the bath degrees of freedom (DOFs).

In this work, we extend a stochastic formulationStockburger and Grabert (2002); Lacroix (2005); Shao (2004) of quantum dissipation by incorporating all three fundamental bath models: non-interacting bosons, fermions and spins. The stochastic Liouville equation (SLE), Eq. (II.4), prescribes a simple yet general form of quantum dissipative dynamics when the bath effects are modelled as colored noises $\xi(t)$ and $\eta(t)$. Different bath models and bath properties are distinguished in the present framework by assigning distinct noise variables and associated statistics. For instance, in dealing with bosonic and fermionic baths, the noises are complex-valued and Grassmann-valued Gaussian processes, respectively, and characterized by the two-time correlation functions such as Eq. (II.4). The Grassmann-valued noises are adopted whenever the environment is composed of fermionic modes as these algebraic entities would bring out the Gaussian characteristics of fermionic modes. For anharmonic environments, such as a spin bath, the required noises are generally non-Gaussian. Two-time statistics cannot fully distinguish these processes and higher order statistics furnished with bath multi-time correlation functions are needed.

Despite the conceptual simplicity of the SLE, achieving stable convergences in stochastic simulations has proven to be challenging in the long-time limit. Even for the most well-studied bosonic bath models, it is still an active research topic to develop efficient stochastic simulation schemsStockburger (2016); Wiedmann et al. (2016); Zhou et al. (2016) today. Our group has successfully applied stochastic path integral simulations to calculate (imaginary-time) thermal distributionsMoix et al. (2012), absorption/emission spectraMoix et al. (2015) and energy transferMoix et al. (2015); Moix and Cao (2013); however, a direct stochastic simulation of real-time dynamics remains formidable. In this study, we consider generic quantum environments that either exhibit non-Gaussian characteristics or invovles ferimonic degrees of freedoms (and associated Grassmann noise in the stochastic formalism). Both scenarios present new challenges to developing efficient stochastic simulations. Hence, in subsequent discussions, all numerical methods developed are strictly deterministic. We note that it is common to derive exact master equationChen and You (2013); Li et al. (2011), hierarchical equation of motionsde Vega (2015); Zhou et al. (2007), and hybrid stochastic-deterministic numerical methodsMoix et al. (2012); Moix and Cao (2013); Zhou et al. (2007) from a stochastic formulation of open quantum theory. In sec. III, we further illustrate the usefulness of our stochastic frmulation by presenting a numerical scheme that would be difficult to obtain within a strictly deterministic framework of open quantum theory. Furthermore, the stochastic formalism gives straightforward prescriptions to compute dynamical quantities such as $\left\langle F_{s}Q_{B}\right\rangle$, which represents system-bath joint observables, as done in a recently proposed theoryYan (2014).

Staring from the SLE, we derive numerical schemes to deterministically simulate quantum dynamics for all three fundamental bath models. The key step is to formally average out the noise variables in the SLE. A common approach is to introduce auxiliary density matrices (ADMs), in close parallel to the hierarchical equation of motion (HEOM) formalismTanimura (2006); *ishizaki_tanimura_jpsp2005, that fold noise-induced fluctuations on reduced density matrix in these auxiliary constructs with their own equations of motions. To facilitate the formal derivation with the noise averaging, we consider two distinct ways to expand the ADMs with respect to a complete set of orthonormal functions in the time domain. In the first case, the basis set corresponds to the eigenfunctions of the bath’s two-time correlation functions. This approach provides an efficient description of open quantum dynamics for bosonic bath models. Unfortunately, it is not convenient to extend this approach to study non-Gaussian bath models. We then investigate another approach inspired by a recent work on the extended HEOM (eHEOM)Tang et al. (2015). In this case, we expand the bath’s multi-time correlation functions in an arbitrary set of orthonormal functions. This approach generalizes eHEOM to the study of non-Gaussian and fermionic bath models with arbitrary spectral densities and temperature regimes. Despite having a slightly more complex form, our fermionic HEOM can be easily related to the existing formalismYe et al. (2016). In this work, we refer to the family of numerical schemes discussed in this work collectively as the generalized hierarchical equations (GHE).

Among the three fundamental bath models, spin baths deserve more attentions. A spin bath can feature very different physical propertiesProkofiev and Stamp, P. C. E. (2000) from the standard heat bath composed of non-interacting bosons; especially, when the bath is composed of localized nuclear / electron spinsHsieh et al. (2012); Kloeffel and Loss (2013), defects and impurities. This kind of spin environment is often of a finite-size and has an extremely narrow bandwidth of frequencies. To more efficiently handle this situation, we consider a dual-fermion mapping that transforms each spin into a pair of coupled fermions. At the expense of introducing an extra set of fermionic DOFs, it becomes possible to recast the non-Gaussian properties of the original spin bath in terms of Gaussian processes in the extended space. In a subsequent work, the paper IIHsieh and Cao (2016), we should further investigate physical properties of spin bath models.

The paper is organized as follows. In Sec. II, we introduce thoroughly the stochastic formalism for open quantum systems embedeed in a generic quantum environment. The SLE is the starting point that we build upon to construct generalized hierarchical equations (GHE), a family of deterministic simulation methods after formally averaging out the noise variables. In Sec. III, we study bosonic baths by expanding the noise processes in terms of the spectral eigenfunctions of the bath’s two-time correlation function. In Sec. IV, we discuss the alternative derivation that generalizes the recently introduced extended HEOM to study non-Gaussian bath models and fermionic bath models. In Sec. V, we introduce the dual-fermion representation and derive an alternative GHE more suitable for spin bath models composed of nuclear / electron spins. A brief summary is given in Sec. VI. In App.A-B, we provide additional materials on the stochastic calculus and Grassmann number to clarify some details of the present work. Appendix C shows how to recover the influence functional theory from the present stochastic formalism. In the last appendix, we present numerical examples to illustrate the methods discussed in the main text.

II A Unified Stochastic Framework For Open Quantum Systems

II.1 The system and bath model

In this study, we consider the following joint system and bath Hamiltonian,

[TABLE]

where the interacting Hamiltonian can be decomposed into factorized forms with $A^{\alpha}$ and $B^{\bar{\alpha}}$ acting on the system and bath, respectively. More specifically, we assume $H_{B}=\sum_{k>0}\omega_{k}b^{{\dagger}}_{k}b_{k}$ and $B^{\alpha}=\sum_{k>0}g_{k}b^{\alpha}_{k}$ where the operator $b^{{\dagger}}_{k}$ can be taken as the bosonic creation operator $a^{{\dagger}}_{k}$, the fermionic creation operator $c^{{\dagger}}_{k}$ or the spin raising operator $\sigma^{+}_{k}$ depending on the specific bath model considered. Throughout the rest of this paper, the system Hamiltonian reads

[TABLE]

where we simply take the system as a spin and the index 0 always refers to the system. Although we adopt a specific system Hamiltonian in Eq. (2), it can be more general.

We consider the factorized initial conditions,

[TABLE]

where two parts are initially uncorrelated and the bath density matrix is commonly taken as a tensor product of thermal states for each individual mode. With Eqs.(1)-(3), the dynamics of the composite system (system and bath) is obtained by solving the von Neumann equation

[TABLE]

If only the system is of interest, then one can trace over the bath DOFs, i.e. $\hat{\rho}_{s}(t)=\text{Tr}_{B}\hat{\rho}(t)$. This straightforward computation (full many-body dynamics then partial trace) soon becomes intractable as the dimension of the Hilbert space scales exponentially with respect to a possibly large number of bath modes.

II.2 Stochastic Decoupling of Many-Body Quantum Dynamics

Many open quantum system techniques have been proposed to avoid a direct computation of Eq. (4). Our starting point is to replace Eq. (4) with a set of coupled stochastic differential equations,

[TABLE]

where the stochastic noises appear as differential Wiener increments, $dW_{\alpha}(t)=\mu_{\alpha}(t)dt$ and $dV_{\alpha}=\nu_{\alpha}(t)dt$ (the white noises $\mu_{\alpha}(t)$ and $\nu_{\alpha}(t)$ will be explicitly defined later). Noting the stochastic relationCao et al. (1996); Lacroix (2005); Diósi et al. (1998) between density matrix and wave function, $\tilde{\rho}_{s}(t)=\left|\psi^{+}(t)\right\rangle\left\langle\psi^{-}(t)\right|$, the two forward/backward wave functions evolves under the time-dependent Hamiltonian, $H^{+}(t)=H_{s}+A\mu(t)$ and $H^{-}(t)=H_{s}+A^{\dagger}\nu^{*}(t)$, respectively. This connection implies an equivalent formulation of stochastic wave function-based theory of open quantum systems with better scalingZhong et al. (2014) for large system simulations. In this work, we focus on the density matrix presentation and attach the tilda symbol on top of all stochastically-evolved density matrices as in Eqs. (5)-(6).

Equation (6) can be further decomposed into corresponding equations for individual modes,

[TABLE]

where $k=1,\dots,N_{B}$, $H_{b,k}=\omega_{k}b^{{\dagger}}_{k}b_{k}$ and $\tilde{\rho}_{B}=\otimes_{k=1}^{N_{B}}\tilde{\rho}_{k}$. The system and bath are decoupled from each other but subjected to the same set of random fields. The stochastic processes ($W_{\alpha}(t)$ and $V_{\alpha}(t)$) in this work can be either complex-valued or Grassmann-valued depending on the bath model under study. To manifest the Gaussian properties of fermionic baths, it is essential to adopt the Grassmann-valued noises. In these cases, it is crucial to maintain the order between fermionic operators and Grassmann-valued noise variables presented in Eqs. (6)-(II.2) as negative signs arise when the order of variables and operators are switched.

The white noises satisfy the standard relations

[TABLE]

where the overlines denote averages over noise realizations. Any other unspecified two-time correlation functions vanish exactly. The order of variables in Eq. (II.2) also matters for Grassmann noises as explained earlier. With these basic properties laid out, we elucidate how to recover Eq. (4) from the stochastic formalism. First, the equation of motion for the joint density matrix $\tilde{\rho}(t)=\tilde{\rho}_{s}(t)\tilde{\rho}_{B}(t)$ is given by

[TABLE]

where the last term is needed to account for all differentials up to $\mathcal{O}(dt)$ as the product of the conjugate pairs of differential Wiener increments such as $\overline{dW_{\alpha}(t)dW^{*}_{\alpha}(t)}=2dt$, contributes a term proportional to $dt$ on average. Taking the noise averages of Eq. (9), the first two terms together yield $-idt[\hat{H}_{s}+\hat{H}_{B},\overline{\tilde{\rho}_{s}\tilde{\rho}_{B}}]$ and the last term gives the system-bath interaction, $-idt[\hat{H}_{I},\overline{\bar{\rho}_{s}\bar{\rho}_{B}}]$. Due to the linearity of the von Neumann equation and the factorized initial condition, the composite system dynamics is given by $\hat{\rho}(t)=\overline{\tilde{\rho}_{s}(t)\tilde{\rho}_{B}(t)}$. To extract the reduced density matrix, we trace out the bath DOFs before taking the noise average,

[TABLE]

In this formulation, it is clear that all the bath-induced dissipative effects are encoded in the trace of the bath’s density matrix.

Because of the non-unitary dynamics implied in the Eqs. (5)-(6), the norm of the stochastically evolved bath density matrices are not conserved along each path of noise realization. The norm conservations only emerge after the noise averaging. This is a common source of numerical instabilities one encounters when directly simulating the simple stochastic dynamics presented so far.

The norm fluctuations of $\tilde{\rho}_{B}(t)$ can be suppressed by modifying the stochastic differential equations above to read,

[TABLE]

where additional stochastic fields

[TABLE]

is introduced to ensure $\text{Tr}_{B}\tilde{\rho}_{B}(t)$ is conserved along each noise path. In Eqs. (11)-(II.2), the top/bottom sign in the symbols $\pm$ ($\mp$) refers to complex-valued/Grassmann-valued noises, respectively. Similarly, $[\cdot,\cdot]_{\mp}$ refers to commutator (complex-valued noise) and anti-commutator (Grassmann-valued noise) in Eq. (11). After these modifications, the exact reduced density matrix of the system is given by $\hat{\rho}_{s}(t)=\overline{\tilde{\rho}_{s}(t)}$.

By introducting $\mathcal{B}^{\alpha}(t)$, we directly incorporate the bath’s response to the random noise in system’s dynamical equation. Equation (11) and the determination of $\mathcal{B}^{\alpha}(t)$ constitute the foundation of open system dynamics in the stochastic framework. In addition to being a methodology of open quantum systems, it should be clear that the present framework allows one to calculate explicitly bath operator involved quantities of interestYan (2014).

II.3 Bath-Induced Stochastic Fields

From Eq. (15), it is clear that $\mathcal{B}^{\alpha}(t)$ can be obtained by formally integrating Eq. (II.2). For simplicity, we take $dV_{\alpha}=dV$ and similarly for $dW_{\alpha}=dW$. This simplification does not compromise the generality of the results presented below and applies to the common case of the spin boson like model when the system-bath interacting Hamiltonian is given by $H_{\text{int}}=\sigma^{z}_{0}\sum_{k>0}g_{k}(b^{+}_{k}+b_{k})$ with $\mathcal{B}(t)=\mathcal{B}^{+}(t)+\mathcal{B}^{-}(t)$.

We first consider the Gaussian baths composed of non-interacting bosons or fermions. The equations of motion for the creation and annihilation operators for individual modes read

[TABLE]

where the top (bottom) sign of $\pm$ should be used when complex-valued (Grassmann-valued) noises are adopted. The expectation values in Eqs. (16)-(17) are taken with respect to $\tilde{\rho}_{k}(t)$. The generalized cumulants are defined as

[TABLE]

For bosonic and fermionic 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/F}(\omega_{k})$ representing either the Bose-Einstein or Fermi-Dirac distribution depending on whether it is a bosonic or fermionic mode.

Replacing the second order cumulants in Eqs. (16)-(17) with an appropriate thermal distribution, one can derive a closed form expression

[TABLE]

where

[TABLE]

stands for the corresponding two-time bath correlation functions. For the Gaussian bath models, Eqs. (11) and (19) together provide an exact account of the reduced system dynamics.

Next we illustrate the treatment of non-Gaussian bath models within the present framework with the spin bath as an example. In the rest of this section, the analysis only applies to the complex-valued noises. The determination of $\mathcal{B}(t)$ still follows the same procedure described at the beginning of this section up to Eqs. (16)-(17). The deviations appear when one tries to compute the time derivatives of the second cumulants. Common to all non-Gaussian bath models, the second cumulants are not time invariant. Instead, by iteratively using Eq. (II.2), the second and higher order cumulants can be shown to obey the following general equation,

[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. The time evolution of these bath cumulants form a simple hierarchical structure with an $n$-th order cumulant influenced directly by the $n+1$-th order cumulants according to the equation above where 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. (18),

[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 do not vanish and persist up to all orders. In any calculations, one should certainly truncate the cumulants at a specific 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. (23). Through this simple prescription, one derives

[TABLE]

where $\Phi_{n,m}(t_{1},\dots,t_{n})$ stands for an $n$-time bath correlation functions with the subscript $m$ used to distinguish the $2^{n-1}$ $n$-time correlation functions appearing in the stochastic integrations (each involves a unique sequence of noise variables) in Eq. (25). Comparing to Eq.(19),one can identify $\Phi_{2,1}(t,s)$ and $\Phi_{2,2}(t,s)$ with $\alpha(t-s)$ and $\alpha^{*}(t-s)$, respectively. Note this derivation assumes the odd-time correlation functions vanish with respect to the initial thermal equilibrium state.

II.4 Stochastic Liouville Equation

At this point, we briefly summarize the unified stochastic formalism. Once $\mathcal{B}(t)$ is fully determined, Eq. (11) can be presented in a simple form, the stochastic Liouville equaiton (SLE),

[TABLE]

where the newly defined color noises are

[TABLE]

In these equations above, the top/bottom signs are associated with complex-valued/Grassmann-valued noises, respectively.

In the cases of bosonic baths, $\mathcal{B}(t)$ is given by Eq. (19) and driven by the complex-valued noises. The color noises then are fully characterized by the statistical properties,

[TABLE]

Several stochastic simulation algorithms have been proposed to solve the SLE with the Gaussian noises.

On the other hand, all previous efforts in the stochastic formulations of fermionic bath end up with derivations of either master equationsZhao et al. (2012); Chen and You (2013) or hierarchicalSuess et al. (2015, 2014) type of coupled equations. The stochastic framework has simply served as a mean to derive deterministic equations for numerical simulations. The lack of direct stochastic algorithm is due to the numerical difficulty to model Grassmann numbers. In this study, we support the view that Grassmann numbers are simply “formal bookkeeping devices” to help formulate the fermionic path integrals and the formal stochastic equations of motion with Gaussian properties. Hence, it is critical to formally eliminate the Grassmann number and the associated stochastic processes, which will be demonstrated in Sec. V with numerical illustrations in App. D.

Going beyond the Gaussian baths, $\mathcal{B}(t)$ is given by Eq. (25) which involves multiple-time stochastic integrals. Formally, one can still use the same definition of noises, Eq. (II.4), and the SLE still prescribes the exact dynamics for the reduced density matrix. The primary factor distinguishing from the Gaussian baths is the statistical characterization of the noises. Higher order statistics are no longer trivial for non-Gaussian processes, and they are determined by the multi-time correlation functions $\Phi_{nm}(t,t_{1},\dots t_{n-1})$ in Eq. (25). For instance, when the fourth order correlation functions are included in the definition of $\mathcal{B}(t)$, additional statistical conditions such as $\overline{\xi(t_{1})\xi(t_{2})\xi(t_{3})\xi(t_{4})}$ would have to be imposed and related to $\{\Phi_{4m}(t_{1},t_{2},t_{3},t_{4}),\Phi_{2m}(t_{i},t_{j})\}$ to fully specify this noise. Since constructing a purely stochastic method to simulate Gaussian processes is already a non-trivial task, simulating non-Gaussian random processes is an even tougher goal.

In the subsequent discussions, we should devise deterministic numerical methods based on the SLE, Eq. (II.4), by formally averaging out the noises. We name the proposed methods in Sec. III - Sec. IV collectively as the generalized hierarchical equations (GHE) in this work. Besides the GHE to be presented, we note sophisticated hybrid algorithmsMoix and Cao (2013); Zhou et al. (2007) could also be constructed to combine advantages of both stochastic and deterministic approaches. We shall leave these potential extensions in a future study.

III Solution I: Spectral Expansion of Stochastic Processes

In this section, we consider a bosonic Gaussian bath. We note Eq. (II.4) encodes the full dissipative effects induced by a bath in the stochastic formalism. The microscopic details, such as Eq. (II.4), of the noise variables become secondary concerns. This observation allows us to substitute any pairs of correlated color noises that satisfy Eq. (II.4) as these statistical conditions alone do not fully specify the noises present in Eq. (II.4). In other words, more than one set of noises can generate identical quantum dissipative dynamics as long as they all satisfy Eq. (II.4) but may differ in other unspecified statistics such as $\overline{\xi(t)\xi^{*}(t^{\prime})}$ and $\overline{\eta(t)\eta^{*}(t^{\prime})}$ etc. This flexibility with the choiceSong et al. (2016) of stochastic processes provides opportunities to fine-tune performances of nuemrical algorithms.

We propose the following decompositionYoung and Grace (2013) of the noise variables

[TABLE]

where ($x_{k}$,$y_{k}$) and ($x^{\prime}_{k}$, $y^{\prime}_{k}$) are independent and real-valued normal variables with mean [math] and variance $1$ while $z_{k}=z^{r}_{k}+iz^{i}_{k}$ are similarly defined but complex-valued. The other unspecified functions are obtained from the spectral expansion of the correlation functions,

[TABLE]

where $\alpha(t-t^{\prime})=\alpha_{R}(t-t^{\prime})+i\alpha_{I}(t-t^{\prime})$, the functions $\chi_{k}(t)$ (complex-valued in general), $\psi_{k}(t)$ (real-valued) and $\phi_{k}(t)$ (real-valued) form independent sets of orthonormal basis function over time domain $[0,T]$ for the simulation. The spectral components can be determined explicitly by solving

[TABLE]

and similar integral equations will yield other sets of basis functions and the associated eigenvalues. The newly defined noises in Eq. (III) can be shown to reproduce the two-time statistics given in Eq. (II.4).

A crucial assumption of the spectral expansion above is that correlation functions should be positive semi-definite. This could be a concern with the quantum correlation functions in the low temperature regime. However, this problem can be addressed by modifying the Hamiltonian and re-define the correlation functions in order to shift the spectral values by a large constant to avoid negative eigenvalues.

We re-label the newly introduced random variables $g_{k}\in\{x_{k},y_{k},x^{\prime}_{k},y^{\prime}_{k},z^{r}_{k},z^{i}_{k}\}$ and expand the reduced density matrix byYoung and Grace (2013)

[TABLE]

where the function $\Phi_{\mathbf{m}}(\mathbf{g})$ with $\mathbf{g}=(g_{1},\dots g_{s})$ is explicitly defined in the second line. The $m$-th Hermite polynomial $\mathcal{H}_{m}(g)$ takes argument of random variables $g$. The total number of random variables $g_{k}$ is given by $s$. Every auxiliary density matrices $\sigma_{\mathbf{m}}(t)$ directly contributes to the determination of $\tilde{\rho}_{s}(t)$.

Substituting Eqs. (III) and (32) into Eq. (II.4) and average over all random variables $g_{k}$, one obtains a set of coupled equations for the density matrices,

[TABLE]

where $\theta_{k}(t)\in\{\sqrt{\lambda_{k}^{R}}\psi_{k}(t),\sqrt{i\lambda^{I}_{k}}\phi_{k}(t),\sqrt{\lambda_{k}/2}\chi_{k}(t),i\sqrt{\lambda_{k}/2}\chi_{k}(t)\}$, the components of $\xi(t)$ in Eq. (III), and similarly $\theta^{\prime}_{k}(t)$ correspond to the components of $\eta(t)$, respectively. In the above equation, $G^{k}_{\mathbf{m}\mathbf{n}}$ is defined by

[TABLE]

where these averages can be done analytically by exploiting the properties of the Hermite polynomials and Gaussian integralsYoung and Grace (2013).

Finally, the exact reduced density matrix is obtained after averaging out random variables in Eq. (32), which can be done by invoking the Gaussian integral identities. The present approach introduces an efficient decomposition of the noise variables and provides an alternative coupling structure for a system of differential equations than the standard HEOM in solving open quantum dynamics. Unfortunately, the present method is not easily generalizable to accommodate the non-Gaussian processes.

IV Solution II: Generalized Hierarchical Equation of Motion

We next present another approach, starting from Eq. (II.4) again, that yields deterministic equations and more easily to accommodate non-Gaussian bath models. This time we utilize Eq. (II.4) as the definitions for the noises, $\xi(t)$ and $\eta(t)$. Following the basic procedure of Ref. Zhou et al., 2007, we average over the noises in Eq. (II.4) to get

[TABLE]

The noise averages yield an auxiliary density matrix (ADM), $\overline{\tilde{\rho}_{s}\xi(t)}=\overline{\rho_{s}\eta(t)}=\overline{\tilde{\rho}_{s}\mathcal{B}(t)}$, by Eq. (II.4). Working out the equation of motion for the ADM, one is then required to define additional ADMs and solve their dynamics too. In this way, a hierarchy of equation of motions for ADMs develops with the general structure

[TABLE]

The time derivatives of $\tilde{\rho}_{s}$ and $\mathcal{B}$ are given by Eqs. (11) and (25), respectively. If we group the ADM’s into a hierarchical tier structure according to the exponent $m$ of $\mathcal{B}^{m}$, then it would be clear soon that the first term of the RHS of Eq. (36) couples the present ADM to ones in the ($m+1$)-th tier, and the last term couples the present ADM to others in the ($m-1$)-th tier.

In the rest of this section, we should materialize these ideas by formulating generalized HEOMs in detail. We separately consider the cases of complex-valued noises (for bosonic and non-Gaussian bath models) and Grassmann-valued noises (for fermionic bath models).

IV.1 Complex-Valued Noise

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

[TABLE]

where $\mathbf{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 cumulant matrices

[TABLE]

where each composed of $2^{n}$ row vectors with indefinite size. 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. (25) can be 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 at the zero-th tier with all $\mathbf{A}_{i}$ being null. Furthermore, the very first ADM we discuss in Eq. (35) can be cast as

[TABLE]

where each ADM 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]

This equation involves 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. (42). We caution that there is no $\mathbf{A}_{0}$ matrix and such a term whenever arises should simply be ignored when interpreting the above equation. After the first term on the RHS of Eq. (IV.1), we only explicitly show the matrices $\mathbf{A}_{n}$ affected in each term of the equation.

This generalized HEOM structure reduces exactly to the recently proposed eHEOMTang et al. (2015) when only $\mathbf{A}_{1}$ cumulant matrix carries non-zero elements, i.e. only the second cumulant expansion of an influence functional is taken into account. It is clear that the higher-order non-linear effects induced by the bath’s $n$-time correlation functions will only appear earliest at the ($n-1$)-th tier expansion.

IV.2 Grassmann-valued Noise

Next we consider the fermionic bath models with Grassmann-valued noises. As discussed earlier, the Grassmann numbers are essential to manifest the Gaussian properties of fermionic baths. In this case, one can significantly simplify the generalized HEOM in the previous section. First, the bath-induced stochastic field can still be decomposed in the form,

[TABLE]

where the indices $m$ and $n$ are suppressed when compared to Eq. (44). This is because the Hamiltonian we consider in this study only allows fermion Gaussian bath model. Each element $a_{j}$ are similarly defined

[TABLE]

where $dU=dW(dV)$ (Grassmann-valued) when $j\leq K/2$ (or $>K/2$). Similar to the algebraic properties of fermionic operators, there are no higher powers of Grassmann numbers and each element $a_{m\mathbf{j}}$ can only appear once. This Pauli exclusion constraint allows us to simplify the representation of fermionic ADMs. We may specify an m-th tier ADM by

[TABLE]

where $\mathbf{n}=(n_{1},\dots n_{K})$ with $n_{i}=0$ or $1$. In this simplified representation, instead of specifying the non-zero elements as in Eq. (45), we layout all elements $a_{mj}$ in an ordered fashion and employ the binary index $n_{i}$ to denote which basis functions contribute to a particular ADM. The tier level of an ADM is determined by the number of basis participating functions, i.e. $m=\sum_{i}n_{i}$.

Following the general procedure outlined in Eq. (36), the m-th tier HEOM reads,

[TABLE]

where $|\mathbf{n}|_{j}=\sum_{i=0}^{j-1}n_{j}$ and $\mathbf{n}_{j,j^{\prime}}$ implies setting $n_{j}=0$ and $n_{j^{\prime}}=1$. In the above equation, $\mathbf{1}_{j}$ is a vector of zero’s except an one at the $j$-th component. The factor such as $n_{j}$ and $(1-n_{j})$ are present to enforce the Pauli exclusion principle associated with the fermions.

The structure of fermionic HEOM certainly resembles that of the bosonic case. However, a few distinctions are worth emphasized. First, it is just the Gaussian bath result including only $\mathbf{A}_{1}$ block matrix when compared to the results in sec. IV.1. Secondly, the fermionic HEOM truncates exactly at some finite number of tiers due to the constraint on the array of binary indices, $n_{i}$. Extra negative signs arise from the permutations to shift the underlying Grassmann-valued stochastic variables from their respective positions in Eq. (IV.2) to the left end of the sequence.

V Spin Bath: Dual-Fermion Transformation

We dedicate an entire section to discuss the spin bath from two distinct perspectives. If one formulates the SLE for a spin bath model in terms of complex-valued noises, then the bath induced stochastic field $\mathcal{B}(t)$ is given by Eq. (25). In this way, the spin bath is a specific example of non-Gaussian bath models. The generalized HEOM formulated earlier can be directly applied in this case. However, in any realistic computations, it is necessary to truncate the statistical characterization of $\mathcal{B}(t)$ up to a finite order of multi-time correlation functions in Eq. (25). While the method is numerically exact, it is computationally prohibitive to calculate beyond the first few higher-order corrections. When the spin bath is large and can be considered as a finite-size approximation to a heat bath, one can show that the linear response approximationMakri (1999) often yields accurate results and a leading order correction should be sufficient whenever needed. The relevance of this leading order correction for spin bath models will be investigated in the paper IIHsieh and Cao (2016).

On the other hand, a spin bath composed of nuclear / electrons spins, as commonly studied in artificial nanostructures at ultralow temperature regime, can beahve very differently from a heat bath composed of non-interacting bosons. There is no particular reason that the linear response and the first leading order correction should sufficiently account for quantum dissipations under all circumstances. In this scenario, it could be useful to map each spin mode onto a pair of coupled fermions. The non-linear mapping allows us to efficiently capture the exact dynamics in an extended Gaussian bath model.

V.1 Dual-Fermion Representation

We consider the following transformation that maps each spin mode into two fermions via,

[TABLE]

where the fermion operators satisfy the canonical anti-commutation relations. One can verify the above mapping reproduces the correct quantum angular momentum commutation relations with the $c$ fermion operators for each spin, while the presence of additional $d$ fermions makes the spin operators associated with different modes commute with each other. We now re-write the Hamiltonian as

[TABLE]

The initial density matrix still maintains a factorized form in the dual-fermion representation,

[TABLE]

where $\hat{\rho}_{c}^{eq}(\beta)$ is the thermal equilibrium state of the $c$ fermions at the original inverse temperature $\beta$ of the spin bath and the $d$ fermions are in the maximally mixed state which is denotes by the identity matrix with dimension $N=2^{n}$ where $n$ is the number of bath modes. A normalization constant is implied to associate with the $\hat{I}_{N}$ matrix.

According to the transformed Hamiltonian in Eq. (53), the system bath coupling now involves three-body interactions, $\sigma^{z}_{0}(d_{k}^{\dagger}+d_{k})(c_{k}^{\dagger}-c_{k})$. Furthermore, the two femrionic baths portrait a non-equilibrium setting with $c$ fermionic bath inherits all physical properties of the original spin bath while $d$ fermionic bath is always initialized in the infinite-temperature limit regardless of the actual state of the spin bath. We first take the system and the d fermions together as an enlarged system and treat the c fermions collectively as a fermionic bath. We introduce the Grassmann noises to stochastically decouple the two subsystems,

[TABLE]

where $B_{k}=g_{k}(c^{\dagger}_{k}-c_{k})$, $\mathcal{B}_{k}=g_{k}\text{Tr}\left\{(c^{\dagger}_{k}+c_{k})\tilde{\rho}_{k}(t)\right\}$, $A_{k}=(d^{\dagger}_{k}+d_{k})$ and $H_{b,k}=\omega_{k}c^{\dagger}_{k}c_{k}$. The density matrices $\tilde{\rho}_{sd}$ denotes the extended system including system spin and all $d$ fermions and $\tilde{\rho}_{k}$ denotes the individual c fermions with $k=1,\dots,N_{B}$. In Eq. (V.1), the noises are defined by

[TABLE]

where $dW_{k}$ and $dV_{k}$ are the standard Grassmann noises defined earlier. The Eq. (55) clearly conserve the norm of $\tilde{\rho}_{k}(t)$ along each noise path, and we will focus on Eq. (V.1) and the stochastic fields, $\mathcal{B}_{k}(t)$.

Our main interest is just the system spin. Hence, we trace out the d fermions in Eq. (V.1) and get

[TABLE]

The auxiliary objects, $\psi^{n}_{\mathbf{k}}$, appearing in Eq. (57) are defined via

[TABLE]

where $k_{1}>\cdots>k_{n}$, (s)Tr either implies standard trace (n is even) or super-trace (n is odd), and the new noises

[TABLE]

A hierarchical structure is implied in Eq. (57), so we derive the equations of motions for the auxiliary objects,

[TABLE]

where $|\mathbf{k}|_{\geq j}\equiv\sum_{k_{i}\geq j}k_{i}$.

Since the spin bath model is mapped onto an effective fermionic problem, the uses of Grassmann noises, Eq. (V.1), will serve as a starting point to develop deterministic numerical methods once the Grassmann noises are integrated out.

V.2 Dual-Fermion GHE

To solve Eq. (V.1), we first define the generalized ADMs

[TABLE]

where $\Omega=(\mathbf{k},\boldsymbol{\alpha})$. The 2 index vectors $\mathbf{k}$ and $\mathbf{j}$ label pairs of coupled $c$ and $d$ fermionic modes; furthermore, the two index vectors are mutually exclusive in the sense a bath mode can appear in just one of the two vectors each time. In this case, the tier-structure of the ADMs are determined by $n+m$. Same as the bosonic and fermionic bath results, the desired reduced density matrix is exactly given by the zero-th tier of ADMs.

Further notational details of Eq. (61) are explained now. The vector $\boldsymbol{\alpha}$ is to be paired with the vector $\mathbf{k}$ to characterize the first set of stochastic fields $\mathcal{B}^{\alpha}_{k}$ in Eq. (61). More precisely, each $(k_{i},\alpha_{i}=\pm)$ labels one of the two possible stochastic fields, $\mathcal{B}^{+}_{k_{i}}=g_{k_{i}}\langle c^{{\dagger}}_{k_{i}}\rangle$ and $\mathcal{B}^{-}_{k_{i}}=g_{k_{i}}\langle c_{k_{i}}\rangle$, associated with $k_{i}$-th c fermion. In dealing with bosonic, fermionic and non-Gaussian baths, we need to explicitly use bath’s multi-time correlation functions via $\mathcal{B}=\sum_{k}\left(\mathbf{B}^{+}_{k_{i}}\pm\mathbf{B}^{+}_{k_{i}}\right)$ when formulating the generalized HEOM approach. In the present case, the stochastic decoupling we introduced in Sec. V.1 dictates that each c fermion acts as a bath and equipped with its own set of stochastic fields as shown in Eq. (V.1). There is no need to expand the bath correlation functions in some orthonormal basis, as each mode’s correlation functions will be treated explicitly in a Fourier decomposition.

Repeat the same steps of the derivation as before, we obtain the generalized HEOM for the spin bath,

[TABLE]

where $\Omega\pm(l,\alpha_{l})$ means a stochastic field $\mathcal{B}^{\alpha_{l}}_{l}$ is either added or removed from the vectors $\mathbf{k}$ and $\boldsymbol{\alpha}$ and, similarly, $\mathbf{j}\pm\mathbf{1}_{l}$ means an index $j=l$ is either added or removed from $\mathbf{j}$.

The range of the index values $\mathbf{k}$ and $\mathbf{j}$ can be extremely large as we explicitly label each microscopic bath modes. Due to the hierarchical structure and the way ADMs are defined, it becomes prohibitively expensive to delve deep down the hierarchical tiers in many realistic calculations. However, the situation might not be as dire as it appears. We already discuss how the present formulation is motivated by the physical spin based environment, such as a collection of nuclear spins in a solid. In such cases, the bath often possess some symmetries allowing simplifications. For instance, most nuclear spins will precess at the same Lamour frequency, and the coupling constant is often distance-dependent. Hence, one can construct spatial “symmetric shells” centered around the system spin in the 3-dimensional real space such that all bath spins inside a shell will more or less share the same frequency and system-bath coupling coefficient. By exploiting this kind of symmetry arguments, one can combine many ADMs defined in Eq. (61) together to significantly reduce the complexity of the hierarchical structures. For a perfectly symmetric bath (i.e. one frequency and one system-bath coupling term), one can use the following compressed ADM,

[TABLE]

where the sum takes into account of all possible combinations of $n$ modes compatible with the requirement that $\sum_{i}\alpha_{i}=s$.

On the other hand, if one deals with a large spin bath described by an effective spectral density then treating the spin bath as an anharmonic environment and usage of the generalize HEOM in Sec. IV.1 will be more appropriate. In fact, in the thermodynamical limit, the spin bath can be accurately approximated as a Gaussian bath and one only needs to invoke $\mathbf{A}_{1}$ block matrix in most calculations.

VI Concluding Remarks

In summary, we advocate the present stochastic framework as a unified approach to extend the study of dissipative quantum dynamics beyond the standard bosonic bath models. We exploit the Itô calculus rule to represent any bilinear interaction between two quantum DOF as white noises. Starting from Eqs. (11)-(II.2), one can derive the SLE, Eq. (II.4), with appropriate statistical conditions, such as Eq. (II.4), that the noises must satisfy. In the Gaussian bath models, the required conditions only involve two-time statistics determined by the bath’s correlation functions. In the case of non-Gaussian bath models, the noises are further characterized by higher order statistics and the multi-time correlation functions.

We devise a family of GHE to solve the SLE with deterministic simulations. We consider two separate orthonormal basis expansions: (1) spectral expansion and (2) generalized HEOM. The spectral expansion, in Sec. III, allows us to solve bosonic bath models efficiently when bath’s two-time correlations assume a simple spectral expansion. This is often the case for correlation functions with a slow decay. The second approach, in Sec. IV, generalizes the eHEOM method to handle multi-time correlation functions in some arbitary set of orthnormal functions. This generalization can provide numerically exact simulations for non-Gaussian (including spin), fermionic and bosonic bath models with arbitray sepctral densities and temperature regimes.

Among the bath models, we extensively discuss the spin bath. When a spin bath is characterized by a well-behaved spectral densityMakri (1999), the generalized HEOM in Sec. IV serves as an efficent approach to simulate dissipative quantum dynamics in a non-Gaussian bath. For situations requiring more than a few higher-order response functions, such as baths composed of almost identical nuclear / electron spins, an alternative approach is to first map the spin bath onto an enlarged Gaussian bath model of fermions via the dual-fermion representation and apply the dual-fermion GHE in Sec. V.2. Numerical examples are illustrated in App. D.

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 Processes

We focus on the complex-valued stochastic processes in this appendix. Additional remarks on Grassmann noises will be made in the following section.

The basic Wiener processes considered in this work is taken to be $W(t)=\int^{t}_{0}ds\nu(s)$ where the complex-valued noise has a mean $\overline{\nu(t)}=0$ and a variance $\overline{\nu(t)\nu^{*}(t^{\prime})}=2\delta(t-t^{\prime})$. Take a uniform discretization of time domain, in each time interval $dt_{i}=t_{i}-t_{i-1}$, each white noise path reduces to a sequence of normal random variables $\{\nu_{1}\dots\nu_{N}\}$. Hence, at each time interval, an identical normal distribution is given,

[TABLE]

The variance is chosen to reproduce the Dirac delta function in the limit $\Delta t\rightarrow 0$. Furthermore, the differential Wiener increments $dW(t)=\nu(t)dt$ satisfy $dW(t)dW^{*}(t)\sim dt$ as required for Brownian motion. The averaging process, implied by the bar on top of stochastic variables, can now be explicitly defined as

[TABLE]

Appendix B Grassmann Numbers and Noises

Grassmann numbers are algebraic constructs that anti-commute among themselves and with any fermionic operators. Given any two Grassmann numbers, $x$ and $y$, and a fermionic operator, $\hat{c}$, they satisfy

[TABLE]

Furthermore, the Grassmann numbers commute with the vacuum state $\left|0\right\rangle$ and, consequently, anti-commute with $\left|1\right\rangle=\hat{c}^{\dagger}\left|0\right\rangle$. Besides the fermionic operators, these numbers commute with everything else such as the bosonic operators and spin Pauli matrices.

Due to the anti-commutativity, there is no higher powers of Grassmann numbers, i.e. $x^{2}=0$. For instance, a single-variate Grassmann function $F(x)=a+bx$ (all variables, a, b, and x, are Grassmann-valued) can only assume this finite Taylor-expanded form. In general every Grassmann function can be decomposed into odd and even parity, $F(x)=F_{\text{even}}(x)+F_{\text{odd}}(x)$. such that

[TABLE]

where $G(x)$ is another arbitrary function with no particular parity assumed. The fermionic thermal equilibrium states are are even-parity Grassmann function when represented in terms of the fermonic coherent states. This even-parity is preserved under linear driving with Grassmann-valued noises. This means all the Grassmann numbers will commute with the fermonic bath density matrices in our study.

Another relevant algebraic property for our study is

[TABLE]

where $x$ is a Grassmann number and sTr$\left\{\cdot\right\}$ is often termed the super-trace.

Finally, we discuss Grassmann-valued white noises. Similar to the discretized complex-valued white noises introduced earlier, we shall take the noise path as a continuum limit of a sequence of Grassmann numbers, $\{x_{i}\}|_{i=1}^{N}$. We will formally treat them as random numbers with respect to Grassmann Gaussians as probability distributions. More precisely, the following integrals yield the desired first two moments (in analogy to the complex-valued normal random variables),

[TABLE]

where the Gaussians should be interpreted by the Taylor expansion: $e^{xx^{\prime}}=1+xx^{\prime}$. In evaluating the integrals above, we recall the standard Grassmann calculus rule that integration with respect to $x$ is equivalent to differentiation with respect to $x$. With these basic set-ups, one can operationally formulate Grassmann noises in close analogy to the complex-valued cases.

Appendix C Relating Stochastic Formalism and Influence Function Theory

The connection between the two formalisms is usually investigated by deriving the stochastic equations from the influence functional theory via the Hubbard-Stratonovich transformationShao (2004); Stockburger and Grabert (2002). Nevertheless, to advocate the stochastic view of quantum dynamics as a rigorous foundation, we establish the connections in the reversed order. We should restrict to the standard bosonic bath models, but extension should be obvious. We first re-write Eq. (5) as

[TABLE]

where the system-bath interaction is given by Eq. (1). In this revised form, it is immediately clear that

[TABLE]

where $T_{\pm}$ is the time-ordering (+) and anti-time-ordering (-) operator. By inserting a complete set of basis $\{\left|\alpha\right\rangle\}$ at each time slice, Eq. (71) can be put in the form,

[TABLE]

where $\boldsymbol{\alpha}=(\alpha,\alpha^{\prime})$ and $\rho(\boldsymbol{\alpha})\equiv\left\langle\alpha^{\prime}\right|\rho\left|\alpha\right\rangle$.

On the other hand, the trace of $\tilde{\rho}_{B}(t)$, governed by Eq. (6), can be expressed as

[TABLE]

where $\mathcal{B}(t)$ is given by Eq. (19). The exact reduced density matrix is then obtained after formally averaging out the noises in the following equation,

[TABLE]

where an explicit evaluation of the noise average on the last line should yield the standard bosonic bath influence functional. To get the influence function, it is useful to contemplate the discretized integrals for the noise averaging,

[TABLE]

where the bath correlation function $C_{i-j}=C(t_{i}-t_{j})$ is given by

[TABLE]

By using the complex-valued Gaussian integral identity,

[TABLE]

the standard Feymann-Vermon influence functional is recovered. The present result is easily generalized when dealing with non-Gaussian baths and $\mathcal{B}(t)$ is potentially characterized by an infinite number of multi-time correlation functions. The noise averaging in this general case will give the cumulant expansion of an influence functional for any bath.

Appendix D Numerical Illustration of Spin Bath Models With Dual-Fermion based HEOM

We present a few numerical results to illustrate the dual-fermion GHE method introduced in this work. Numerical examples with generalized HEOM approach will be further studied in a separate work, the paper IIHsieh and Cao (2016). We will consider various cases of a pure dephasing model,

[TABLE]

An analytical expression for the off-diagonal matrix element of the reduced density matrix reads,

[TABLE]

with

[TABLE]

where $\Omega_{k}=\omega_{k}\sqrt{1+(2g_{k}/\omega_{k})^{2}}$.

First, we consider a 50-spin bath with the parameters $(\omega_{k},g_{k})$ sampled from the discretization of an Ohmic bath. We use the general dual-fermion GHE scheme, Eq. (V.2), to simulate the off-diagonal matrix element for the density matrix. Figure 1 shows the results in the weak coupling (panel a) and the strong coupling (panel b) cases. Due to each spin is modelled as a bath, it becomes prohibitive to delve into further tiers. Nevertheless, with a shallow 2-tier hierarchy, the results seem to do reasonably well in the short-time limit.

In the second case, we consider the spin star modelSinayskiy et al. (2012); *wang_guo_eurphys13; *breuer_burgarth_prb04; *petruccione_prb06 where all the bath spins look identical, i.e. $\omega_{k}=\omega^{\prime}$ and $g_{k}=g^{\prime}$. This is an often used model to analyze spin bath models. As shown by the results in Fig. 2, it is critical to go deep down the hierarchical tiers in order to recover the correct quantum dissipations. One can only generate this many tiers through compressing the auxiliary density matrices as in Eq. (63). This second example illustrates the kind of scenarios where dual-fermion GHE could provide an accurate account of quantum dynamics induced by a spin bath.