# Spin relaxation dynamics with a continuous spin environment: the dissipaton
equation of motion approach

Wenxiang Ying<EMAIL_ADDRESS>Department of Chemistry, University of
Rochester, 120 Trustee Road, Rochester, New York 14627, USA  Yu Su  Department
of Chemical Physics, University of Science and Technology of China, Hefei,
Anhui 230026, China  Zi-Hao Chen  Department of Chemical Physics, University
of Science and Technology of China, Hefei, Anhui 230026, China  Yao Wang
<EMAIL_ADDRESS>Department of Chemical Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  Pengfei Huo
<EMAIL_ADDRESS>The Institute of Optics, Hajim School of
Engineering, University of Rochester, Rochester, New York 14627, USA
Department of Chemistry, University of Rochester, 120 Trustee Road, Rochester,
New York 14627, USA

###### Abstract

We present the quantum dynamics of a spin coupling to a bath of independent
spins via the dissipaton equation of motion (DEOM) approach. The bath,
characterized by a continuous spectral density function, is composed of spins
that are independent level systems described by the $\mathfrak{su}(2)$ Lie
algebra. This represents an extreme class of anharmonic environment. Based on
the conclusion drawn by Suarez and Silbey [J. Chem. Phys. 95, 9115 (1991)] and
Makri [J. Chem. Phys. 111, 6164 (1999)] that the spin bath can be mapped to a
Gaussian environment under its linear response limit, we derive the
fluctuation-dissipation theorem (FDT) of the spin bath from a microscopic
perspective, and generalize the discussion to the case of arbitrary bath spin
quantum number $S$. Next, the time-domain Prony fitting decomposition scheme
is applied to the bare-bath time correlation function (TCF) given by FDT to
generate the exponential decay basis (or pseudo modes) for DEOM construction.
The accuracy and efficiency of this strategy has been justified by a variety
of numerical results. We envision this work provides new insights to extend
the hierarchical equations of motion (HEOM) and DEOM approach to certain types
of anharmonic enviroments with arbitrary TCF or spectral density.

##  I Introduction

The dynamics of a two-level-system (TLS) coupling to a dissipative environment
has been extensively investigated over the past decades. The most widely and
systematically studied one is the spin-boson system,[1, 2, 3, 4] which has
various applications in physics and chemical dynamics in condensed phase[1, 2,
3, 4, 5, 6]. The environmental part of the spin-boson system is a set of non-
interacting harmonic oscillators. This is reasonable as Caldeira and Leggett
[7] had justified the universality of bosonic heat baths consisting of an
infinite number of harmonic oscillators which are linearly coupled to the
system in the real physical world.

Apart from the bosonic environment, another typical environment of interest is
a bath consisting of a set of spins, [8, 9] which can be regarded as an
extreme example of anharmonic environment. The spin-spin-bath (SSB) model
denotes a TLS coupled with a dissipative spin bath. It has recently drawn
tremendous attention due to the gradually accumulating phenomena in physical
setups. Here are some evidences. It is predicted and justified that at very
low temperature, the dynamics of magnetic nanomolecules, such as
$\mathrm{Fe}_{8}$, $\mathrm{Mn}_{12}$, is strongly influenced by nuclear
spins; [10, 11, 12, 13, 14] for solid-state quantum computing devices whose
qubits are typically electron spins, such as GaAs quantum dot [15, 16, 17] and
diamonds with nitrogen-vacancy (NV) centres [18, 19], it is inevitable for the
qubits to be coupled to environmental spins. The spin environment is also
concerned in stylized quantum measurement setups, [20, 21] the studies of
quantum phase transition, [17, 22, 23] and more recently, the radical pair
spin relaxation as well as its applications in quantum information processing,
[24, 25, 26, 27] etc. Anharmonic but very simple, the SSB model provides new
insights for quantum dynamics in condensed phase.

The theoretical treatment to the anharmonic spin bath has first been discussed
by Suarez and Silbey [28], Makri [29], and later by Yan [30]. They proved that
under the thermodynamics limit (also known as the linear response limit), the
generally anharmonic spin bath approximately satisfies the Gaussian
statistics. As a result, the spin bath can be effectively characterized by the
familiar boson bath with a modified spectral density function, which is much
easier to be treated by a wide range of quantum dynamics methods. On this
basis, higher order nonlinear effects can be further studied by considering a
finite number of bath spins. Previous work on the quantum dynamics of SSB
types of models have been carried out extensively by iterative path integral
approach based on influence functional, [29, 31, 32] the multilayer
multiconfiguration time-dependent Hartree (ML-MCTDH) [33] and the closely
related surrogate Hamiltonian [34, 35] approach, polaron-transformed master
equation, [36] Nakajima-Zwanzig type of generalized master equation,[37] and
the generalized hierarchical equations of motion (gHEOM) that based on the
stochastic Liouville equation with perturbative expansion scheme, [38, 39]
etc.

A variety of approximation-based methods are also developed and applied for
the spin bath problems, such as the time-dependent perturbation approach,[40]
phase space quasi-classical methods, [41] etc. A number of cluster expansion
(CE) approaches [42], including the cluster correlation expansion [43, 44,
45], linked cluster expansion [46], and the associated dynamical mean field
theory [47] have been developed to study the many-body bath time evolution.
Nakamura and Tanimura [48] studied the dynamics of a TLS that interacts with a
subenvironment consisting of a one-dimensional XXZ spin chain using the
hierarchical Schrödinger equations of motion (HSEOM), despite that the noise
generated from the spin lattice is non-Gaussian and non-local. [2, 49] The
series of work done by Fay, Lindoy, and Manolopoulos [50, 51, 24, 25, 52, 53]
had provided a systematic theoretical framework for electron spin relaxation
in radical pairs that exposed to an environment of nuclear spins, which
includes full quantum mechanical treatment of all spin degrees of freedom
using tensor network propagation strategy, and master equation approaches
based on the Schulten-Wolynes semiclassical treatment to the nuclear spins.

Despite the fruitful progresses, these approaches all have their specific
limitations. For example, the direct path integral-based method are
numerically efficient only for short memory length; [54, 55, 56] the
numerically exact ML-MCTDH approach relies on a discretization strategy for
the continuous bath spectral density, is computational costly to reach
numerical convergence for models with high bath cutoff frequency, so that
Born-Oppenheimer (B-O) type of approximation needs to be made; [1, 33] MCTDH
is even more computational costly under the finite-temperature case due to its
dependence on the Monte Carlo sampling strategy. [57] Besides, almost all of
the work mentioned above only take care of the spin-1/2 case. In realistic
situations, however, the nuclear spin quantum numbers are usually much larger
than $1/2$.

Regarding to the limitations of current studies on the spin bath that
mentioned above, a more general, numerically efficient and accurate
theoretical framework needs to be developed. The dissipaton equation of motion
(DEOM) is a statistical quasi-particle theory for quantum dissipative
dynamics. Featured by the very powerful dissipaton algebras, it not only just
recovers HEOM formalism, [58, 59, 60] but also identifies the physical
meanings of the dynamical variables. [9, 61] Equipped with the time-domain
Prony fitting decomposition ($t$-PFD) [62] scheme, by which the environmental
TCF can be accurately decomposed into exponential sums, the efficiency and
applicability of HEOM/DEOM is significantly improved. It is notable that
$t$-PFD scheme is, in principle, applicable to arbitrary bath TCF or spectral
density. Under the linear response limit, the theory of effective spectral
density function provides a perfect platform for us to play with $t$-PFD.

In this work, by taking advantage of the linear response limit, we generalize
the theory of effective spectral density function for spin bath with arbitrary
bath spin quantum number $S$. Then we apply DEOM to study the spin relaxation
dynamics of a TLS coupling to the effective bosonic environment. The bath TCF
decomposition is achieved by $t$-PFD. This paper is organized as follows. In
Section II, we briefly review the spin bath model, its linear response
properties, bosonic DEOM, and $t$-PFD scheme; in Section III, we present
various numerical testing results for SSB models, including zero- and finite-
temperature, weak and strong coupling, unbiased and biased cases with
comparisons against ML-MCTDH; in Section IV, we briefly summarized the major
advantage of our method.

##  II Model and methodology

The total TLS-plus-bath composite Hamiltonian of an open quantum system reads
as

$\displaystyle\hat{H}=\hat{H}_{\mathrm{S}}+\hat{h}_{\textsc{B}}+\hat{H}_{\mathrm{SB}},$
(1)

where
$\hat{H}_{\mathrm{S}}\equiv\epsilon\hat{\sigma}_{z}+\Delta\hat{\sigma}_{x}$ is
the system Hamiltonian of a TLS with energy bias $\epsilon$ and off-diagonal
coupling $\Delta$,
$\hat{\bm{\sigma}}\equiv(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})$
denotes to the Pauli matrices; $\hat{h}_{\textsc{B}}$ is the bath Hamiltonian
which usually consists of a macroscopic number of noninteracting particles
that can be bosons, fermions and/or spins. $\hat{H}_{\mathrm{SB}}$ carries the
system-bath interaction. The influence of the bath entails quantum statistical
mechanics description, in which the thermodynamics limit is naturally assumed.
We set $\hbar\equiv 1$ throughout this paper for the sake of convenience.

###  II.1 The spin bath models

The spin bath and its interaction with the system can be described by the
Hamiltonian below, [29, 33, 39]

$\displaystyle\hat{h}_{\textsc{B}}$
$\displaystyle=\sum_{j=1}^{N}\omega_{j}\hat{s}^{j}_{z},$ (2a)
$\displaystyle\hat{H}_{\mathrm{SB}}$
$\displaystyle=\sum_{a}\hat{Q}_{a}\otimes\hat{F}_{a}.$ (2b)

The bath Hamiltonian $\hat{h}_{\textsc{B}}$ describes $N$ independent spins
which are distinguishable. It is diagonal in the
$\\{(\hat{\bm{s}}^{j})^{2},\hat{s}^{j}_{z}\ |\ j=1,\cdots,N\\}$ eigen
representation, where $\hat{\bm{s}}^{j}$ are the spin operators associated
with the $j$-th bath mode, characterized by energy difference $\omega_{j}$.
They form the $\otimes_{j=1}^{N}\mathfrak{su}(2)$ Lie algebra,
$[\hat{s}^{i}_{\alpha},\hat{s}^{j}_{\beta}]=i\delta_{ij}\epsilon_{\alpha\beta\gamma}\hat{s}^{i}_{\gamma}$,
where $\alpha,\beta,\gamma$ denotes to the Cartesian components of the spin
matrices, $\epsilon_{\alpha\beta\gamma}$ is the 3-D Levi-Civita tensor,
$\delta_{ij}$ is the Kronecker symbol. $\\{\hat{Q}_{a}\\}$ and
$\\{\hat{F}_{a}\\}$ denotes to the system and bath dissipation modes,
respectively.

There are extensive types of spin-spin interaction, such as the Ising type,
[63] the Heisenberg type, [64] etc., whose general expression takes the form
$\hat{\sigma}_{\alpha}\hat{s}^{j}_{\beta}$. Among the various choices, the
most common types can be summarized as [39]

$\displaystyle\hat{H}_{\mathrm{SB}}=\begin{cases}\frac{1}{\sqrt{S}}\sum_{j=1}^{N}c_{j}\hat{s}^{j}_{x}\hat{\sigma}_{z},\vspace{0.15cm}\\\
\frac{1}{\sqrt{S}}\sum_{j=1}^{N}c_{j}\hat{s}^{j}_{z}\hat{\sigma}_{z},\vspace{0.15cm}\\\
\frac{1}{\sqrt{S}}\sum_{j=1}^{N}c_{j}(\hat{s}^{j}_{x}\hat{\sigma}_{x}+\hat{s}^{j}_{y}\hat{\sigma}_{y}),\vspace{0.15cm}\\\
\frac{1}{\sqrt{S}}\sum_{j=1}^{N}c_{j}(\hat{s}^{j}_{x}\hat{\sigma}_{x}+\hat{s}^{j}_{y}\hat{\sigma}_{y}+\hat{s}^{j}_{z}\hat{\sigma}_{z}),\end{cases}$
(3)

where $c_{j}$ are the coupling coefficients between the system dissipation
operators and the $j$-th bath mode, $S$ is the spin quantum number of the bath
spins. In this paper, we shall mainly focus on the first interaction form in
Eq.(3). According to Caldeira and Leggett, [7] the bath as well as its
coupling to the system might be described by the spectral density function,
defined as

$\displaystyle J_{ab}(\omega)$
$\displaystyle\equiv\frac{\pi}{2}\sum_{j=1}^{N}c^{*}_{aj}c_{bj}\delta(\omega-\omega_{j}),$
(4)

where the coupling coefficients are assumed to obey the general scaling rule:
$c_{aj}\sim 1/\sqrt{N}$.

It is worth noting that the spin bath is not generally a Gaussian environment.
Suarez and Silbey [28], Makri [29] have remarkably shown that under the limit
of $N\rightarrow+\infty$, the SSB model with the first interaction form in
Eq.(3) and bath spin quantum number $S=1/2$, can be rigorously mapped onto the
familiar spin-boson model with an effective spectral density,

$\displaystyle[J_{ab}]_{\mathrm{eff}}(\omega;\beta)=J_{ab}(\omega)\tanh(\beta\omega/2),$
(5)

with the subscript $a=b$ for single mode case. More specifically, only the
second order term remains in the cumulant expansion of the influence
functional as $N\rightarrow+\infty$. The higher order cumulants
$\mathcal{O}(1/\sqrt{N})$ will disappear under the scaling limit. As a result,
the noise spectrum of the bath dissipation operators are Gaussians, [28]
leading to the linear response limit. 111There is a straightforward example
with respect to the pure dephasing case, in which $\Delta=0$ such that the
dephasing dynamics has analytical solution, as discussed by Rao and Kurizki
[95], Hsieh and Cao [39]

###  II.2 Linear response of the spin bath, and the generalized theory of
effective spectral density function

As is discussed above that the continuous spin bath approximately satisfies
Gaussian statistics. The influence of Gaussian environments are completely
characterized with the linear response functions of hybrid bath modes in the
isolated bare-bath subspace. For bosonic environment, it can be defined via
the commutator as

$\displaystyle\chi_{ab}(t)\equiv
i\langle[\hat{F}_{a}(t),\hat{F}_{b}(0)]\rangle_{\textsc{B}}.$ (6)

Here
$\hat{F}(t)=e^{i\hat{h}_{\textsc{B}}t}\hat{F}(0)e^{-i\hat{h}_{\textsc{B}}t}$
exerts the stochastic force, and
$\langle(\cdot)\rangle_{\textsc{B}}\equiv\mathrm{tr}_{\textsc{B}}[(\cdot)\rho^{eq}_{\textsc{B}}(T)]$
denotes ensemble average in the bath subspace, with
$\rho^{eq}_{\textsc{B}}(T)\equiv
e^{-\beta\hat{h}_{\textsc{B}}}/\mathrm{tr}_{\textsc{B}}[e^{-\beta\hat{h}_{\textsc{B}}}]$.
For fermionic environments, the similar concept can also be defined via the
anti-commutator, [2, 9]
$G_{ab}(t)\equiv\langle\\{\hat{F}_{a}(t),\hat{F}_{b}(0)\\}\rangle_{\textsc{B}}$,
which is known as the single-particle Green’s function. The causality Fourier
transform of $\chi_{ab}(t)$ and $G_{ab}(t)$ is defined as

$\displaystyle\chi_{ab}(\omega)\equiv\int_{0}^{\infty}dt\ e^{i\omega
t}i\langle[\hat{F}_{a}(t),\hat{F}_{b}(0)]\rangle_{\textsc{B}},$ (7a)
$\displaystyle G_{ab}(\omega)\equiv\int_{0}^{\infty}dt\ e^{i\omega
t}\langle\\{\hat{F}_{a}(t),\hat{F}_{b}(0)\\}\rangle_{\textsc{B}}.$ (7b)

One can easily check their symmetry. The spectral density functions can be
evaluated using time-reversal symmetry (TRS) as [2]

$\displaystyle J_{ab}(\omega)$
$\displaystyle\equiv\mathrm{Im}[\chi_{ab}(\omega)]=\frac{1}{2}\int_{-\infty}^{+\infty}dt\
e^{i\omega t}\langle[\hat{F}_{a}(t),\hat{F}_{b}(0)]\rangle_{\textsc{B}},$ (8a)
$\displaystyle J^{\prime}_{ab}(\omega)$
$\displaystyle\equiv\mathrm{Re}[G_{ab}(\omega)]=\frac{1}{2}\int_{-\infty}^{+\infty}dt\
e^{i\omega t}\langle\\{\hat{F}_{a}(t),\hat{F}_{b}(0)\\}\rangle_{\textsc{B}},$
(8b)

corresponding to the bosonic and fermionic cases, respectively. Likewise, one
can define the TCF as well as its spectrum functions as

$\displaystyle C_{ab}(t)$
$\displaystyle\equiv\langle\hat{F}_{a}(t)\hat{F}_{b}(0)\rangle_{\textsc{B}},$
(9a) $\displaystyle C_{ab}(\omega)$ $\displaystyle\equiv\int_{0}^{\infty}dt\
e^{i\omega t}\langle\hat{F}_{a}(t)\hat{F}_{b}(0)\rangle_{\textsc{B}}.$ (9b)

Without loss of generality, we consider the case of linear system-bath
coupling with only one dissipation mode,

$\displaystyle\hat{F}\equiv\frac{1}{\sqrt{S}}\sum_{j}c_{j}\hat{s}^{j}_{x},$
(10)

for the spin bath with $S=1/2$. Here we drop all the subscripts for
simplicity. One will finally arrive at its FDT with respect to the auto-TCF of
the bath dissipation operator,

$\displaystyle
C(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}d\omega\frac{e^{-i\omega
t}J^{\prime}(\omega)}{1+e^{-\beta\omega}}.$ (11)

See Appendix A for detailed derivations. On the other hand, one can
independently obtain that:

$\displaystyle C(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}d\omega\
\frac{e^{-i\omega t}J_{\mathrm{eff}}(\omega;\beta)}{1-e^{-\beta\omega}},$ (12)

where $J_{\mathrm{eff}}(\omega;\beta)\equiv
J(\omega)=J^{\prime}(\omega)\tanh(\beta\omega/2)$, recovering Eq.(5), the
well-known result. So that very interestingly, the continuous spin environment
is isomorphic to a boson environment with temperature-dependent effective
spectral density function, whose zero-temperature limit gives rise to
$J^{\prime}(\omega)=J(\omega)$. For this reason, the zero-temperature spin-
boson model is also widely known as taking the spin-bath limit. 222 Note that
historically, Schotte[96] first pointed out that at low temperature and long
time limit ($t\gg\omega_{c}^{-1}$), the dynamics given by the Kondo
Hamiltonian in a bosonic picture is the same as that of the corresponding
fermionic operators. This pattern also emerges in the theoretical validation
of the Caldeira-Leggett model [7, 8]

Due to the generality of the linear response limit, the theory of effective
spectral density is widely applicable, so that it is not restricted to $S=1/2$
case, but with arbitrary $S$, as long as the system-bath coupling is linear.
We also provide discussions upon arbitrary spin $S$ case in the Appendix A,
which can be viewed as a generalization of the effective spectral density
function theory.

###  II.3 Bosonic DEOM formalism

Based on the previous discussion, the open quantum system problem with a
continuous spin environment will be exactly mapped to the familiar spin-boson
type of model with an effective spectral density function, which has the
Hamiltonian description as below,

$\displaystyle\hat{H}=\hat{H}_{\mathrm{S}}+\frac{1}{2}\sum_{j}\omega_{j}(\hat{x}^{2}_{j}+\hat{p}^{2}_{j})+\sum_{a}\hat{Q}_{a}\otimes\sum_{j}c^{\prime}_{aj}\hat{x}_{aj},$
(13)

where $\hat{x}_{j}$, $\hat{p}_{j}$ are the conjugated coordinate-momentum
pairs that satisfy the Heisenberg commutation relations, $\\{\hat{Q}_{a}\\}$
are the original system dissipation mode, and $\\{c^{\prime}_{aj}\\}$ are the
rescaled coupling coefficients due to the effective spectral density. The
problem is now readily to be solved by HEOM/DEOM.

The DEOM theory is a statistical quasi-particle theory for quantum dissipative
dynamics, describing the influence of bulk environments using only a few
number of quasi-particles, the dissipatons. [9, 61] They arise strictly from
the linear bath coupling component:

$\displaystyle\hat{F}_{a}=\sum_{k=1}^{K}\hat{f}_{ak},$ (14)

with single-damping parameters given by

$\displaystyle\langle\hat{f}_{ak}(t)\hat{f}_{bj}(0)\rangle_{\textsc{B}}=\delta_{kj}\eta_{abk}e^{-\gamma_{ak}t},$
(15a)
$\displaystyle\langle\hat{f}_{bj}(0)\hat{f}_{ak}(t)\rangle_{\textsc{B}}=\delta_{kj}\eta_{ab\overline{k}}^{*}e^{-\gamma_{ak}t}.$
(15b)

The associated index $\overline{k}$ in Eq.(15b) is defined as
$\gamma_{a\overline{k}}=\gamma_{ak}^{*}$ to preserve TRS. Further denote

$\displaystyle\langle\hat{f}_{ak}\hat{f}_{bj}\rangle^{>}_{\textsc{B}}\equiv\langle\hat{f}_{ak}(0^{+})\hat{f}_{bj}(0)\rangle_{\textsc{B}}=\delta_{kj}\eta_{abk},$
(16a)
$\displaystyle\langle\hat{f}_{bj}\hat{f}_{ak}\rangle^{<}_{\textsc{B}}\equiv\langle\hat{f}_{bj}(0)\hat{f}_{ak}(0^{+})\rangle_{\textsc{B}}=\delta_{kj}\eta_{ab\overline{k}}^{*},$
(16b)

for later use in the dissipaton algebra. Note that they are different from
$\langle\hat{f}_{ak}\hat{f}_{bj}\rangle_{\textsc{B}}$. Eq.(14) and (15) leads
to

$\displaystyle\langle\hat{F}_{a}(t)\hat{F}_{b}(0)\rangle_{\textsc{B}}=\sum_{k=1}^{K}\eta_{abk}e^{-\gamma_{ak}t},$
(17)

and its complex conjugation.

The dynamical variables in DEOM are the dissipaton density operators (DDOs):

$\displaystyle\rho^{(n)}_{\mathbf{n}}(t)\equiv\mathrm{Tr}_{\textsc{B}}\left[\Big{(}\prod_{ak}\hat{f}^{n_{ak}}_{ak}\Big{)}^{\circ}\rho_{\mathrm{T}}(t)\right],$
(18)

where $\rho_{\mathrm{T}}(t)$ is the time-dependent total density matrix, the
product of dissipatons inside $(\cdots)^{\circ}$ means irreducible. And
$(\mathrm{c\text{-}number})^{\circ}=0$. Bosonic dissipatons are symmetric
under permutation,
$(\hat{f}_{ak}\hat{f}_{bj})^{\circ}=(\hat{f}_{bj}\hat{f}_{ak})^{\circ}$. Each
DDO in Eq.(18) represents a specific configuration of
$\mathbf{n}\equiv\\{\cdots,n_{ak},\cdots|\ a=1,\cdots,M;\ k=1,\cdots,K\\}$,
with $n=\sum_{ak}n_{ak}$ dissipatons in total (i.e., the number of tiers). We
also denote that the associated DDO’s index $\bm{\mathrm{n}}^{\pm}_{ak}$
differs from $\bm{\mathrm{n}}$ at the specified $n_{ak}$ by $\pm 1$, which
means $n_{ak}$ is replaced by $n_{ak}\pm 1$.

The DEOM formalism can be constructed according to the dissipaton algebra,
which includes the generalized diffusion equation and generalized Wick’s
theorem. The generalized diffusion equation arises from the single-damping
character in Eq.(15), that

$\displaystyle\mathrm{Tr}_{\textsc{B}}\left[\left(\frac{\partial\hat{f}_{ak}}{\partial
t}\right)_{\textsc{B}}\rho_{\mathrm{T}}(t)\right]=-\gamma_{ak}\mathrm{Tr}_{\textsc{B}}\left[\hat{f}_{ak}\rho_{\mathrm{T}}(t)\right].$
(19)

The generalized diffusion equation is applicable for the
$\hat{h}_{\textsc{B}}$-action,

$\displaystyle\rho^{(n)}_{\bm{\mathrm{n}}}(t;h^{\times}_{\textsc{B}})$
$\displaystyle\equiv\mathrm{Tr}_{\textsc{B}}\left\\{\Big{(}\prod_{ak}\hat{f}^{n_{ak}}_{ak}\Big{)}^{\circ}[\hat{h}_{\textsc{B}},\rho_{\mathrm{T}}(t)]\right\\}$
$\displaystyle=\mathrm{Tr}_{\textsc{B}}\left\\{\Big{[}\Big{(}\prod_{ak}\hat{f}^{n_{ak}}_{ak}\Big{)}^{\circ},\hat{h}_{\textsc{B}}\Big{]}\
\rho_{\mathrm{T}}(t)\right\\}$
$\displaystyle=-i\left(\sum_{ak}n_{ak}\gamma_{ak}\right)\rho^{(n)}_{\bm{\mathrm{n}}}(t),$
(20)

where $h^{\times}_{\textsc{B}}\ \cdot\equiv[\hat{h}_{\textsc{B}},\ \cdot]$,
the second line of Eq.(II.3) arises from the equivalence between the
Schrödinger and Heisenberg prescription, and the last line goes with
Heisenberg equations of motion, $(\partial\hat{f}_{ak}/\partial
t)_{\textsc{B}}=-i[\hat{f}_{ak},\hat{h}_{\textsc{B}}]$. Eq.(II.3) summarizes
the contribution by the bath Hamiltonian to the DDOs dynamics.[9, 61]

The generalized Wick’s theorem deals with the system-hybrid-bath interaction,
reading as

$\displaystyle\mathrm{Tr}_{\textsc{B}}\left[\Big{(}\prod_{ak}\hat{f}^{n_{ak}}_{ak}\Big{)}^{\circ}\hat{f}_{bj}\rho_{\mathrm{T}}(t)\right]$
$\displaystyle=\rho^{(n+1)}_{\bm{\mathrm{n}}^{+}_{bj}}(t)$ (21a)
$\displaystyle\ \ \
+\sum_{ak}n_{ak}\langle\hat{f}_{ak}\hat{f}_{bj}\rangle^{>}_{\textsc{B}}\rho^{(n-1)}_{\bm{\mathrm{n}}^{-}_{ak}}(t),$
$\displaystyle\mathrm{Tr}_{\textsc{B}}\left[\Big{(}\prod_{ak}\hat{f}^{n_{ak}}_{ak}\Big{)}^{\circ}\rho_{\mathrm{T}}(t)\hat{f}_{bj}\right]$
$\displaystyle=\rho^{(n+1)}_{\bm{\mathrm{n}}^{+}_{bj}}(t)$ (21b)
$\displaystyle\ \ \
+\sum_{ak}n_{ak}\langle\hat{f}_{ak}\hat{f}_{bj}\rangle^{<}_{\textsc{B}}\rho^{(n-1)}_{\bm{\mathrm{n}}^{-}_{ak}}(t).$

They will be used in evaluating the commutator action of linear system-bath
coupling terms. The bosonic DEOM formalism is now readily to be constructed,
reading as [9, 61]

$\displaystyle\dot{\rho}^{(n)}_{\bm{\mathrm{n}}}(t)$
$\displaystyle=-\left(i\mathcal{L}_{\mathrm{S}}+\sum_{ak}n_{ak}\gamma_{ak}\right)\rho^{(n)}_{\bm{\mathrm{n}}}(t)-i\sum_{ak}\mathcal{Q}^{\times}_{a}\rho^{(n+1)}_{\bm{\mathrm{n}}^{+}_{ak}}(t)$
$\displaystyle\ \ \
-i\sum_{abk}n_{ak}\left(\eta^{\prime}_{abk}\mathcal{Q}^{\times}_{b}+i\eta^{{}^{\prime\prime}}_{abk}\mathcal{Q}^{\circ}_{b}\right)\rho^{(n-1)}_{\bm{\mathrm{n}}^{-}_{ak}}(t).$
(22)

The involved superoperators and coefficients are defined as below:

$\displaystyle\mathcal{L}_{\mathrm{S}}\hat{O}\equiv[\hat{H}_{\mathrm{S}},\hat{O}],\
\ \ \ \ \mathcal{Q}^{\times}_{a}\hat{O}\equiv[\hat{Q}_{a},\hat{O}],\ \ \ \ \
\mathcal{Q}^{\circ}_{a}\hat{O}\equiv\\{\hat{Q}_{a},\hat{O}\\},$
$\displaystyle\ \ \ \ \ \ \ \ \
\eta^{\prime}_{abk}\equiv\frac{\eta_{abk}+\eta^{*}_{ab\overline{k}}}{2},\ \ \
\ \
\eta^{{}^{\prime\prime}}_{abk}\equiv\frac{\eta_{abk}-\eta^{*}_{ab\overline{k}}}{2i}.$

The RK-4/RK-45 algorithm [67] are usually adopted as the numerical propagation
scheme of Eq.(II.3). There are also on-the-fly filtering algorithms available
for acceleration. [68]

###  II.4 Time-domain Prony fitting decomposition

In this context, the central problem in DEOM is to decomposite the bare-bath
TCF into a sum of exponential series. Based on FDT, this can be realized by
expanding-over-pole strategies, such as Matsubara spectral decomposition (MSD)
[2] and Pad$\mathrm{\acute{e}}$ spectral decomposition (PSD) [69, 70, 71], or
various least-square fitting schemes [72, 73, 74, 75]. The traditional
expanding-over-poles strategies are usually restricted to certain forms of
bath spectral density. Here we choose the $t$-PFD strategy. The resulting
numerical efficiency of HEOM/DEOM is optimized to a great extent, especially
in low temperature regimes that are usually inaccessible for other methods.

The $t$-PFD scheme is intrinsically based on least-squares fitting algorithm,
in which the real and imaginary parts of TCF are fitted separately. It is easy
to extract from the bosonic FDT (c.f. Eq.(12)) that

$\displaystyle\mathrm{Re}[C(t)]$
$\displaystyle=\frac{1}{\pi}\int_{0}^{+\infty}d\omega\
J_{\mathrm{eff}}(\omega;\beta)\coth\left(\frac{\beta\omega}{2}\right)\cos(\omega
t)$ $\displaystyle=\frac{1}{\pi}\int_{0}^{+\infty}d\omega\
J(\omega)\cos(\omega t),$ (23a) $\displaystyle\mathrm{Im}[C(t)]$
$\displaystyle=\frac{1}{\pi}\int_{0}^{+\infty}d\omega\
J_{\mathrm{eff}}(\omega;\beta)\sin(\omega t)$
$\displaystyle=\frac{1}{\pi}\int_{0}^{+\infty}d\omega\
J(\omega)\tanh(\beta\omega/2)\sin(\omega t),$ (23b)

where the second line in Eq.(23a) and (23b) are the explicit results for
$S=1/2$ case. Intriguingly, temperature dependence is only carried by the
imaginary part of TCF, in line with fermion bath but the reverse of boson
bath.

Next, we target at optimizing $K=K_{r}+K_{i}$ in

$\displaystyle\mathrm{Re}[C(t)]=\sum_{k=1}^{K_{r}}\zeta_{k}e^{-\lambda_{k}t},\
\ \
\mathrm{Im}[C(t)]=\sum_{k=1}^{K_{i}}\zeta^{\prime}_{k}e^{-\lambda^{\prime}_{k}t},$
(24)

where $K_{r}$ and $K_{i}$ are the number of terms in real and imaginary part
fitting, respectively. We accordingly denote the $t$-PFD strategy as
$K_{r}+K_{i}$. See Ref. 62 and reference therein for detailed procedures of
$t$-PFD, as well as the numerical benchmarks for several commonly used
spectral density functions. We also provide more illustrations and examples on
the performances of TCF fitting using $t$-PFD in Appendix B, with respect to
the most challenging SSB models we tested.

##  III Computational Details, Results and discussions

In this section, we present the numerical benchmark results of the bosonic
DEOM equipped with $t$-PFD for various SSB models. The total Hamiltonian is as
below:

$\displaystyle\hat{H}_{\mathrm{S}}$
$\displaystyle=\epsilon\hat{\sigma}_{z}+\Delta\hat{\sigma}_{x},$ (25a)
$\displaystyle\hat{h}_{\textsc{B}}$
$\displaystyle=\sum_{j=1}^{N}\omega_{j}\hat{s}^{j}_{z},$ (25b)
$\displaystyle\hat{H}_{\mathrm{SB}}$
$\displaystyle=\hat{\sigma}_{z}\otimes\sum_{j}\sqrt{2}c_{j}\hat{s}^{j}_{x},$
(25c)

We use the Ohmic form [8] with exponential cutoff for the description of the
continuous bath as well as its interaction with the system,

$\displaystyle J(\omega)=\frac{\pi}{2}\alpha\omega e^{-\omega/\omega_{c}}.$
(26)

In this expression, $\alpha$ is the Kondo parameter that characterizes the
system-bath coupling strength, and $\omega_{c}$ is the bath cut-off frequency.
Since Eq.(26) is defined via Eq.(4), it is proved to be in line with Eq.(8b),
according to the discussions conducted in II.2 and Appendix A.

The DEOM propagation uses the fourth-order Runge-Kutta (RK-4) integrator with
time step of 0.0025/$\Delta$, together with the on-the-fly filtering algorithm
[68] with given error tolerance for acceleration.

###  III.1 Zero-temperature spin relaxation dynamics and localization

Figure 1: Population dynamics of the zero-temperature SSB models. The spin
baths are parameterized as (a) $\alpha=0.5$, $\omega_{c}/\Delta=1$. (b)
$\alpha=0.1$, $\omega_{c}/\Delta=6$. (c) $\alpha=0.2$, $\omega_{c}/\Delta=10$.
(d) $\alpha=0.5$, $\omega_{c}/\Delta=10$. (e) $\alpha=0.75$,
$\omega_{c}/\Delta=10$. (f) $\alpha=0.5$, $\omega_{c}/\Delta=40$. The ML-MCTDH
results are digitized from Ref. 33.

Fig. 1 presents the population dynamics
$P(t)\equiv\langle\hat{\sigma}_{z}(t)\rangle$ for zero-temperature unbiased
SSB models (with $\epsilon=0$, to keep in consistence with the previous work).
Even though they equal to the spin-boson models, it provides us with the first
glimpse on the power of $t$-PFD. Here we uniformly set the error tolerance of
on-the-fly filtering algorithm as $5\times 10^{-7}$. We have also confirmed
that the “empirical standard” error tolerance of $1\times 10^{-5}$ remains
accurate enough in the present cases when $\alpha\leq 0.5$. For the model of
Fig. 1e with coupling strength $\alpha=0.75$, the convergence test suggests
that the error tolerance should be set no larger than $1\times 10^{-6}$. All
the numerically converged results are obtained with $t$-PFD strategy $5+5$
(Fig. 1a - d, f) or $6+5$ (Fig. 1e) to ensure the accuracy in TCF fitting, and
large enough number of tiers in the hierarchic expansion (here 20 will be
satisfying). The population dynamics obtained by DEOM (black solid lines) are
compared to ML-MCTDH (red dots). In all of the models presented, DEOM with
$t$-PFD perfectly match the ML-MCTDH results, including the Rabi oscillations
and the incoherent relaxations.

Another important phenomenon about the unbiased SSB/spin-boson model at zero
temperature is the localization, [8] that the population dynamics quickly
reaches a biased stationary value. It is a typical phenomenon when the time
scale of the bath is comparable to or longer than that of the subsystem. [76,
77, 78, 79] Within the adiabatic or intermediate between adiabatic and
nonadiabatic regime (that with a modest $\omega_{c}/\Delta$ value), a large
coupling strength ($\alpha>1$) will bring about a large barrier height along
the adiabatic double-well potential energy surface, such that localization of
the population can be induced. [33]

Fig. 2 presents the convergence test of the localization model (with
$\omega_{c}/\Delta=1,\ \alpha=10$) using different $t$-PFD strategies. This
computation is rather challenging, as it generally requires a lot of memory
and CPU time to reach convergence. To ensure the numerical accuracy, we turn
off on-the-fly filtering module, and the number of tiers is set as $45$ to
stay accurate enough. As is observed that $t$-PFD strategy $2+2$ is already
able to capture the localization phenomenon, but not accurate enough for
population dynamics. At the expense of greater computational cost, the result
generated with $t$-PFD strategy $3+3$ and $4+4$ exhibits better accuracy. It
takes more than $110$ hours of CPU time (Intel Xeon Gold $6330$ CPU @$2.00$
GHz with $36$ cores) with memory requirement no less than $300$ GB to produce
the $4+4$ curve. The convergence is in good agreement with ML-MCTDH. See also
Appendix B for details about the accuracy of TCF fitting under different
$t$-PFD strategies.

Figure 2: Convergence test of the localization model ($\omega_{c}/\Delta=1,\
\alpha=10$) using $t$-PFD strategies $2+2$, $3+3$ and $4+4$. Comparisons are
made against ML-MCTDH (digitized from Ref. 33).

###  III.2 Finite-temperature spin relaxation dynamics, and the localization-
delocalization phase transition

Figure 3: Temperature-dependence in the coherent-incoherent transition of
population dynamics. The spin bath models are parametrized as (a)
$\omega_{c}/\Delta=6,\ \alpha=0.5$. (b) $\omega_{c}/\Delta=10,\ \alpha=0.75$.
(c) $\omega_{c}/\Delta=1,\ \alpha=10$. The DEOM results (solid lines) are
compared to the ML-MCTDH results (solid dots, digitized from Ref. 33).

Next, we turn to the finite-temperature cases. Fig. 3 presents the
temperature-dependence of spin relaxation dynamics for several representative
yet challenging SSB models.

In Fig. 3a, $P(t)$ under different temperatures, with bath parameters
$\omega_{c}/\Delta=6,\ \alpha=0.5$ are shown (corresponding to Fig. 5a of Ref.
33). All the curves are computed after reaching a good convergence using
$t$-PFD strategy $5+5$. Fig. 3b is similar to Fig. 3a but more challenging,
with $\omega_{c}/\Delta=10,\ \alpha=0.75$ (corresponding to Fig. 6c of Ref.
33). So we take $t$-PFD strategy $6+5$ to reach the numerical convergence,
even though $5+5$ will already be accurate enough. The comparisons are also
made against ML-MCTDH. One can observe that at zero-temperature, the
population dynamics exhibits clearly an incoherent decay to the equilibrium
value $P=0$. Increasing temperature by a little bit might greatly change the
paradigm of decay; further increasing temperature will induce stronger
coherent motion in the short time period. The temperature susceptibility is
also closely related to the coupling strength $\alpha$. [33] It is well-known
that at zero-temperature, the coherent-incoherent boundary is at $\alpha=0.5$
in the scaling limit. [8] For the SSB model, high temperature can
significantly raise up this boundary value. As is also reported by Shao and
Hänggi [80], Forsythe and Makri [81]. It is very counter-intuitive yet
interesting that higher temperature will slightly enhance the coherence for
the central TLS in the nonadiabatic regime ($\omega_{c}/\Delta\gg 1$).
Physical interpretation of this abnormal phenomenon can be made from the
analysis of the mapping spin-boson model with effective spectral density:
$J_{\mathrm{eff}}(\omega)$ has a smaller magnitude when temperature increases,
which renders weaker coupling to the central spin and win over the thermal
excitation quenching effect, resulting in more coherent dynamics for the
system TLS. [33] Our results are in good agreement with the previous work.

Fig. 3c presents the temperature-dependence of the localization model
($\omega_{c}=1,\ \alpha=10$). All the curves are computed using $t$-PFD
strategy $4+4$ with good convergence. As is observed that when temperature
increases, the poplation distribution gets a faster decay and the bias
gradually disappears, which is actually a kind of phase transition. [8, 82]
This is because increasing temperature will decrease the relaxation time of
the bath, so that the system TLS become delocalized.

Although our results show DEOM with $t$-PFD are overall in excellent agreement
with the ML-MCTDH results, we shall point out that for a certain number of
finite-temperature models presented here, especially for the one with
$k_{B}T/\Delta=0.1$ and $0.167$ in Fig. 3b, there still exist minor
discrepancies between DEOM and ML-MCTDH results even by eye inspection. This
should be reasonable because ML-MCTDH adopts discretized bath modes as well as
B-O type of approximation for high frequency bath modes, [33] it is also
susceptible to the tensor network propagation scheme; while HEOM/DEOM is able
to use rigorously continuous bath modes. The major resource of error for
HEOM/DEOM comes from the accuracy of bath TCF fitting after reaching
convergence. Another possible reason lies in the finite-temperature strategy.
ML-MCTDH adopts Monte Carlo importance sampling techniques [83] to evaluate
the thermal Boltzmann operator, which is usually hard to reach numerical
convergence; while HEOM/DEOM resorts to FDT, being a deterministic pathway.

###  III.3 More results on the finite-temperature biased models

As we know that in realistic situations, the system is not always unbiased.
For example, when a Zeeman field is applied to the central spin, the energy
degeneracy will be broken, which is a common experimental set up to study the
radical pair spin relaxation dynamics.[24] To this reason, we further study
the spin relaxation dynamics for biased models under finite-temperature, and
compare them to the corresponding spin-boson model with same bath parameters.

Figure 4: Population dynamics and von Neumann entropy of the biased SSB and
spin-boson models. The boson/spin baths are parameterized as (a) $\alpha=0.4$,
$\omega_{c}/\Delta=1$, $\beta\Delta=0.25$. (b) $\alpha=0.4$,
$\omega_{c}/\Delta=2$, $\beta\Delta=1.0$. (c) $\alpha=0.4$,
$\omega_{c}/\Delta=2$, $\beta\Delta=5.0$.

Fig. 4 shows the converged population dynamics as well as von Neumann entropy
for several biased TLS (with $\epsilon=\Delta$) interacting with a spin bath
(SSB) or a boson bath (spin-boson). The von Neumann entropy reads as

$\displaystyle
S_{\mathrm{vN}}(t)\equiv-\mathrm{Tr}[\hat{\rho}_{\mathrm{S}}(t)\ln\hat{\rho}_{\mathrm{S}}(t)],$
(27)

where $\hat{\rho}_{\mathrm{S}}(t)$ is the density matrix of the system TLS at
time $t$. Three different models from high temperature ($\beta=0.25$) to low
temperature ($\beta=5$) are presented here. As is observed that for high
temperature model in Fig. 4a, the SSB model keeps much better quantum
coherence than the corresponding spin-boson model. The von Neumann entropy
growth is also slower. As expected, when temperature decreases, the behavior
of the SSB model gradually agrees with the spin-boson model.

##  IV Concluding Remarks

We present the numerical benchmark results of the spin relaxation dynamics of
various SSB models with continuous spectral density function by using DEOM
with $t$-PFD. The idea of mapping boson bath with temperature-dependent
effective spectral density is validated from a microscopic perspective, and
generalized to arbitrary bath spin quantum number $S$. The highlight of this
strategy is that one can approximate a class of non-Gaussian bath to the
Gaussian one under the linear response limit. By applying the $t$-PFD to the
bath TCF, one obtains the exponential decay basis to construct DEOM and
propagate. The accuracy and numerical efficiency of this strategy is
illustrated by various of examples. DEOM with $t$-PFD provides an excellent
agreement with the ML-MCTDH results reported in literature, despite the rather
challenging model parameters. Moreover, comparing to previous research work on
the SSB model based on ML-MCTDH, modular path integral, etc., DEOM equipped
with $t$-PFD has the advantage that the bath spectral density is rigorously
continuous, no B-O type of approximation is needed; and the finite-temperautre
strategy is based on FDT rather than stochastic sampling, thus allowing long
time propagation. In summary, it serves as a novel and more efficient
benchmark scheme for quantum dynamics with a spin bath.

The presented theoretical work provide new insights to extend HEOM/DEOM to
certain types of non-Gaussian environments with arbitrary bath TCFs or
spectral density functions. Future research work on methods development shall
be carried out on more complicated types of interacting system-bath models
that listed in Eq.(3), with incorporation of quadratic or higher order system-
bath couplings, as well as the nonlinear response effects; we would also like
to seek for applications in the simulation of radical pair spin relaxation
dynamics that exposed to an enviroment of nuclear spins.

###### Acknowledgements.

This work was supported by the National Science Foundation CAREER Award under
Grant No. CHE-1845747. Computing resources were provided by the Center for
Integrated Research Computing (CIRC) at the University of Rochester. W. Y.
would like to sincerely appreciate YiJing Yan for his meticulous guidance when
he was at USTC, and for his valuable comments on this manuscript.

## Data Availability

The data that support the findings of this work are available from the
corresponding author under reasonable request.

##  Appendix A FDT for the spin bath, and the effective bosonic environment

###  A.1 $S=1/2$ case

For spin-1/2 particles, we adopt the well-known Jordan-Wigner correspondence
[84] for the spin operators:

$\displaystyle\hat{d}_{j}^{\dagger}$
$\displaystyle\equiv\hat{s}^{j}_{x}+i\hat{s}^{j}_{y},$ (28a)
$\displaystyle\hat{d}_{j}$
$\displaystyle\equiv\hat{s}^{j}_{x}-i\hat{s}^{j}_{y},$ (28b)
$\displaystyle\hat{s}^{j}_{z}$
$\displaystyle=\hat{d}_{j}^{\dagger}\hat{d}_{j}-\frac{1}{2}.$ (28c)

where $\hat{d}_{j}$, $\hat{d}_{j}^{\dagger}$ are the mapping fermionic
annihilation/creation operators with anti-commutation relation on the same
site. However, on different sites we have bosonic commutation relations, which
means spins on different sites commute, unlike fermions which anti-commute.
For this reason, spins are also referred to as hard-core bosons or spinless
fermions, [85] whose commutation relations can be summarized as

$\displaystyle\\{\hat{d}_{j},\hat{d}_{j}\\}=\\{\hat{d}_{j}^{\dagger},\hat{d}_{j}^{\dagger}\\}=0,\
\ \ \\{\hat{d}_{j},\hat{d}_{j}^{\dagger}\\}=1;$ (29a)
$\displaystyle\\{\hat{d}_{j},\hat{d}_{k}\\}=2\hat{d}_{j}\hat{d}_{k},\ \ \
\\{\hat{d}^{\dagger}_{j},\hat{d}^{\dagger}_{k}\\}=2\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k},$
$\displaystyle\\{\hat{d}_{j},\hat{d}^{\dagger}_{k}\\}=2\hat{d}_{j}\hat{d}^{\dagger}_{k},\
\ \ k\neq j;$ (29b)
$\displaystyle[\hat{d}_{j},\hat{d}_{j}]=[\hat{d}_{j}^{\dagger},\hat{d}_{j}^{\dagger}]=0,\
\ \ [\hat{d}_{j},\hat{d}_{j}^{\dagger}]=1-2\hat{d}_{j}^{\dagger}\hat{d}_{j};$
(29c)
$\displaystyle[\hat{d}_{j},\hat{d}_{k}]=[\hat{d}_{j}^{\dagger},\hat{d}_{k}^{\dagger}]=[\hat{d}_{j},\hat{d}_{k}^{\dagger}]=0,\
\ \ k\neq j.$ (29d)

As is shown that the hard-core bosons possess characters of both bosons and
fermions, one can either define the linear response functions using commutator
or the single particle Green’s function using anti-commutator. Here we will
show both possibilities. For convenience, let’s first write down the time
evolution of single bath operators,

$\displaystyle\hat{d}_{j}^{\dagger}(t)$
$\displaystyle=e^{i\hat{h}_{\textsc{B}}t}\hat{d}_{j}^{\dagger}(0)e^{-i\hat{h}_{\textsc{B}}t}=\hat{d}_{j}^{\dagger}(0)e^{i\omega_{j}t},$
$\displaystyle\hat{d}_{j}(t)$
$\displaystyle=e^{i\hat{h}_{\textsc{B}}t}\hat{d}_{j}(0)e^{-i\hat{h}_{\textsc{B}}t}=\hat{d}_{j}(0)e^{-i\omega_{j}t},$
(30)

where we have used the well-known Baker-Campbell-Hausdorff identity. [86]

We will start with the anti-commutator version. For each pair of mapping
spinless fermionic creation/annihilation operators $\hat{d}_{j}^{\dagger}$ and
$\hat{d}_{j}$ that satisfy the commutation relations defined in Eq.(29), the
single particle Green’s functions can be defined and evaluated as

$\displaystyle\langle\\{\hat{d}_{j}(t),\hat{d}^{\dagger}_{k}(0)\\}\rangle_{\textsc{B}}$
$\displaystyle=\delta_{jk}e^{-i\omega_{j}t},$
$\displaystyle\langle\\{\hat{d}^{\dagger}_{j}(t),\hat{d}_{k}(0)\\}\rangle_{\textsc{B}}$
$\displaystyle=\delta_{jk}e^{i\omega_{j}t},$
$\displaystyle\langle\\{\hat{d}_{j}(t),\hat{d}_{k}(0)\\}\rangle_{\textsc{B}}$
$\displaystyle=\langle\\{\hat{d}^{\dagger}_{j}(t),\hat{d}^{\dagger}_{k}(0)\\}\rangle_{\textsc{B}}=0.$
(31)

So we have

$\displaystyle\langle\\{\hat{F}(t),\hat{F}(0)\\}\rangle_{\textsc{B}}=\sum_{j}c^{2}_{j}\cos(\omega_{j}t).$
(32)

As a result, the spectral density function is evaluated as (c.f. Eq.(8b))

$\displaystyle J^{\prime}(\omega)$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{+\infty}dt\ e^{i\omega
t}\langle\\{\hat{F}(t),\hat{F}(0)\\}\rangle_{\textsc{B}}$
$\displaystyle=\frac{\pi}{2}\sum_{j}c^{2}_{j}[\delta(\omega-\omega_{j})+\delta(\omega+\omega_{j})],$
(33)

which is an extension of the result given by Caldeira and Leggett [7] (see
also Eq.(4)) to negative frequencies $\omega<0$, while ensuring it to be an
even function.

The derivation via the commutator version is very similar to the anti-
commutator one. Following the same procedure, one can easily obtain the bare-
bath linear response function as

$\displaystyle
i\langle[\hat{F}(t),\hat{F}(0)]\rangle_{\textsc{B}}=\sum_{j}c^{2}_{j}\langle(1-2\hat{d}_{j}^{\dagger}\hat{d}_{j})\rangle_{\textsc{B}}\sin(\omega_{j}t).$
(34)

For independent spin $S=1/2$ particles, we have

$\displaystyle
Z_{\textsc{B}}\equiv\mathrm{Tr}_{\textsc{B}}[e^{-\beta\hat{h}_{\textsc{B}}}]=\prod_{j}2\cosh\left(\frac{\beta\omega_{j}}{2}\right),$
(35)

one immediately obtains

$\displaystyle
i\langle[\hat{F}(t),\hat{F}(0)]\rangle_{\textsc{B}}=\sum_{j}c^{2}_{j}\tanh\left(\frac{\beta\omega_{j}}{2}\right)\sin(\omega_{j}t).$
(36)

And the corresponding spectral density function reads as (c.f. Eq.(8a))

$\displaystyle J(\omega)$
$\displaystyle\equiv\frac{1}{2}\int_{-\infty}^{+\infty}dt\ e^{i\omega
t}\langle[\hat{F}(t),\hat{F}(0)]\rangle_{\textsc{B}}$ (37)
$\displaystyle=\frac{\pi}{2}\sum_{j}c^{2}_{j}\tanh\left(\frac{\beta\omega_{j}}{2}\right)[\delta(\omega-\omega_{j})-\delta(\omega+\omega_{j})],$

which is an odd function. On the other hand, we can rewrite Eq.(37) as

$\displaystyle J(\omega)$
$\displaystyle=\frac{\pi}{2}\sum_{j}c^{2}_{j}\tanh\left(\frac{\beta\omega}{2}\right)[\delta(\omega-\omega_{j})+\delta(\omega+\omega_{j})]$
$\displaystyle=J^{\prime}(\omega)\tanh\left(\frac{\beta\omega}{2}\right),$
(38)

where $J^{\prime}(\omega)$ is defined in Eq.(A.1), giving rise to the
temperature-dependent effective spectral density.

Combining Eq.(32) and (36), one obtains

$\displaystyle C(t)$
$\displaystyle\equiv\langle\hat{F}(t)\hat{F}(0)\rangle_{\textsc{B}}=\frac{1}{2}\sum_{j}c^{2}_{j}\left[\frac{e^{-i\omega_{j}t}}{1+e^{-\beta\omega_{j}}}+\frac{e^{i\omega_{j}t}}{1+e^{\beta\omega_{j}}}\right]$
$\displaystyle=\frac{1}{2}\sum_{j}c^{2}_{j}[\delta(\omega-\omega_{j})+\delta(\omega+\omega_{j})]\frac{e^{-i\omega
t}}{1+e^{-\beta\omega}}$
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{+\infty}d\omega\frac{e^{-i\omega
t}J^{\prime}(\omega)}{1+e^{-\beta\omega}},$ (39)

which is just Eq.(11) in the main text, the FDT for spin bath. It has an
equivalent bosonic FDT formalism if one takes
$J_{\mathrm{eff}}(\omega;\beta)\equiv
J(\omega)=J^{\prime}(\omega)\tanh(\beta\omega/2)$, giving rise to Eq. (12).
Thus, the spin-boson problem with an effective spectral density function
arises naturally.

An alternative but similar argument can be done by using the coupled-fermion
representation for spin operators ($S=1/2$), discussed by Mattis, et al.,[87,
88, 39] reading as

$\displaystyle\hat{s}^{j}_{+}$
$\displaystyle=\hat{c}_{j}^{\dagger}(\hat{d}_{j}+\hat{d}_{j}^{\dagger}),$
(40a) $\displaystyle\hat{s}^{j}_{-}$
$\displaystyle=(\hat{d}_{j}+\hat{d}_{j}^{\dagger})\hat{c}_{j},$ (40b)
$\displaystyle\hat{s}^{j}_{z}$
$\displaystyle=\hat{c}_{j}^{\dagger}\hat{c}_{j}-\frac{1}{2},$ (40c)

where $\hat{c},\hat{c}^{\dagger}$ and $\hat{d},\hat{d}^{\dagger}$ are two sets
of fermion operators that anti-commute with each other. On the other hand,
since the mapping fermionic bath operators still commute with the system
operators rather than anti-commute, it should also lead to the same bosonic
DEOM formalism, as is studied by Jin, et al. [89] based on the influence
functional.

Figure 5: $t$-PFD fitting results for the spin bath with
$\omega_{c}/\Delta=1,\ \alpha=10$. (a), (b) real and imaginary part fitting
results for the zero-temperature model using 2, 3, 4 and 5 terms,
respectively. (c), (d) real and imaginary part fitting results for different
finite-temperature models using 4 terms; the dashed lines represent the exact
TCF, and the solid lines represent the fitting results, which are almost
overlapped.

###  A.2 Generalization to arbitrary spin $S$

Adopting the same angular momentum raising/lowering operators that are defined
as

$\displaystyle\hat{s}^{j}_{+}\equiv\hat{s}^{j}_{x}+i\hat{s}^{j}_{y},\ \ \ \ \
\hat{s}^{j}_{-}\equiv\hat{s}^{j}_{x}-i\hat{s}^{j}_{y},$ (41)

so that the original $\otimes_{j=1}^{N}\mathfrak{su}(2)$ Lie algebra becomes

$\displaystyle[\hat{s}^{i}_{+},\hat{s}^{j}_{-}]=2\hat{s}^{i}_{z}\delta_{ij},\
\ \ \ \ [\hat{s}^{i}_{z},\hat{s}^{j}_{\pm}]=\pm\hat{s}^{i}_{\pm}\delta_{ij}.$
(42)

Their time dependence can be evaluated as

$\displaystyle\hat{s}^{j}_{\pm}(t)=e^{i\hat{h}_{\textsc{B}}t}\hat{s}^{j}_{\pm}(0)e^{-i\hat{h}_{\textsc{B}}t}=\hat{s}^{j}_{\pm}(0)e^{\pm
i\omega_{j}t}.$ (43)

The bare-bath partition function can be evaluated as

$\displaystyle
Z_{\textsc{B}}=\prod_{j}\left(2\sum_{k=0}^{[S]}\cosh\left((k+\\{S\\})\beta\omega_{j}\right)-\Delta(S)\right),$
(44)

where $[S]$, $\\{S\\}$ are the integer and fractional part of $S$, with
$S=[S]+\\{S\\}$; $\Delta(S)=1$ for $S$ being integers, and $\Delta(S)=0$ for
$S$ being half integers. As a result,

$\displaystyle\langle\hat{s}^{j}_{z}\rangle_{\textsc{B}}=-\frac{\partial\ln
Z_{\textsc{B}}}{\partial(\beta\omega_{j})},\ \ \ \ \
\langle(\hat{s}^{j}_{z})^{2}\rangle_{\textsc{B}}=\frac{\partial^{2}\ln
Z_{\textsc{B}}}{\partial(\beta\omega_{j})^{2}}.$ (45)

It is also easy to obtain the bare-bath single particle Green’s function and
linear response function as:

$\displaystyle\langle\\{\hat{F}(t),\hat{F}(0)\\}\rangle_{\textsc{B}}$
$\displaystyle=\frac{1}{S}\sum_{j}c^{2}_{j}[S(S+1)-\langle(\hat{s}^{j}_{z})^{2}\rangle_{\textsc{B}}]\cos(\omega_{j}t),$
(46a) $\displaystyle i\langle[\hat{F}(t),\hat{F}(0)]\rangle_{\textsc{B}}$
$\displaystyle=-\frac{1}{S}\sum_{j}c^{2}_{j}\langle\hat{s}^{j}_{z}\rangle_{\textsc{B}}\sin(\omega_{j}t).$
(46b)

And the corresponding spectral density functions can be evaluated by Eq.(8a)
and (8b) as

$\displaystyle J^{\prime}(\omega)$
$\displaystyle=\frac{\pi}{2S}\sum_{j}c^{2}_{j}[S(S+1)-\langle(\hat{s}^{j}_{z})^{2}\rangle_{\textsc{B}}]$
$\displaystyle\ \ \ \ \
\times[\delta(\omega-\omega_{j})+\delta(\omega+\omega_{j})],$ (47a)
$\displaystyle J(\omega)$
$\displaystyle=-\frac{\pi}{2S}\sum_{j}c^{2}_{j}\langle\hat{s}^{j}_{z}\rangle_{\textsc{B}}[\delta(\omega-\omega_{j})-\delta(\omega+\omega_{j})].$
(47b)

They are derived from a microscopic perspective, satisfying the correct
symmetry; however, neither of them give rise to the original definition of
bath spectral density function that given by Caldeira and Leggett [7]. One
might choose to still using Eq.(4) and directly getting $J(\omega)$ in
Eq.(47b) as the effective spectral density, or do modifications to make it in
accordance with Eq.(47a) or (47b). For example, if we use (c.f. Eq.(47a)):

$\displaystyle J^{\prime}(\omega)$
$\displaystyle\equiv\frac{\pi}{2S}\sum_{j}c^{2}_{j}[S(S+1)-\langle(\hat{s}^{j}_{z})^{2}\rangle_{\textsc{B}}]\Big{|}_{\omega_{j}=\omega}$
$\displaystyle\ \ \ \ \
\times[\delta(\omega-\omega_{j})+\delta(\omega+\omega_{j})],$ (48)

where the footnote $\omega_{j}=\omega$ means replacing all the $\omega_{j}$ in
the prefactors by $\omega$. Then the theory of effective spectral density
function can arise as

$\displaystyle C(t)$
$\displaystyle\equiv\langle\hat{F}(t)\hat{F}(0)\rangle_{\textsc{B}}=\frac{1}{\pi}\int_{-\infty}^{+\infty}d\omega
e^{-i\omega t}\frac{J_{\mathrm{eff}}(\omega;\beta,S)}{1-e^{-\beta\omega}},$
$\displaystyle J_{\mathrm{eff}}(\omega;\beta,S)\equiv
J(\omega)=J^{\prime}(\omega)\zeta(\omega;\beta,S),$ (49)

with

$\displaystyle\zeta(\omega;\beta,S)$
$\displaystyle\equiv\frac{1-e^{-\beta\omega}}{2}$ (50) $\displaystyle\ \ \
\times\frac{S(S+1)-\langle(\hat{s}^{j}_{z})^{2}\rangle_{\textsc{B}}-\langle\hat{s}^{j}_{z}\rangle_{\textsc{B}}}{S(S+1)-\langle(\hat{s}^{j}_{z})^{2}\rangle_{\textsc{B}}}\Bigg{|}_{\omega_{j}=\omega}.$

One can check that $\zeta(\omega;\beta,S=1/2)=\tanh(\beta\omega/2)$.

On the other hand, under the high spin limit of $S\gg 1$, Eq.(47b) will reduce
to the conventional form,

$\displaystyle
J(\omega)=\frac{\pi}{2}\sum_{j}c^{2}_{j}[\delta(\omega-\omega_{j})-\delta(\omega+\omega_{j})],$
(51)

which is in line with the original definition in Eq.(4) but with odd
analytical continuation. As a result, the spin bath under the high spin limit
is isomorphic to the boson bath. This can be understood via the Holstein-
Primakoff transformation, [90]

$\displaystyle\hat{s}^{i}_{+}$ $\displaystyle=\sqrt{2S}\
\sqrt{1-\frac{\hat{b}^{\dagger}_{i}\hat{b}_{i}}{2S}}\
\hat{b}_{i}\approx\sqrt{2S}\ \hat{b}_{i},$ $\displaystyle\hat{s}^{i}_{-}$
$\displaystyle=\sqrt{2S}\ \hat{b}^{\dagger}_{i}\
\sqrt{1-\frac{\hat{b}^{\dagger}_{i}\hat{b}_{i}}{2S}}\approx\sqrt{2S}\
\hat{b}^{\dagger}_{i},\ $ $\displaystyle\hat{s}^{i}_{z}$
$\displaystyle=S-\hat{b}^{\dagger}_{i}\hat{b}_{i},$ (52)

where $\hat{b}^{\dagger}_{i}$ and $\hat{b}_{i}$ are bosonic
creation/annihilation operators that satisfy the Heisenberg commutation
relations. Consequently,

$\displaystyle\hat{F}=\frac{1}{\sqrt{S}}\sum_{j}c_{j}\hat{s}^{j}_{x}\approx\sum_{j}\frac{c_{j}}{\sqrt{2}}(\hat{b}_{j}+\hat{b}^{\dagger}_{j})\equiv\sum_{j}c_{j}\hat{x}_{j},$
(53)

where $\hat{x}_{j}\equiv(\hat{b}_{j}+\hat{b}^{\dagger}_{j})/\sqrt{2}$. Eq.(53)
recovers the bath dissipation operator of the conventional spin-boson model.

One should find it straightforward to generalize the discussions above to more
complicated interacting system-bath models, which might contain multiple
dissipation modes as is listed in Eq.(3), by carrying out the very similar
arguments. The spectral density functions will be anistropic in such cases.
Even more general, the linear response limit can be easily applied to general
finite baths [91, 28], with multiple dissipation modes. i.e., the individual
bath particles are general level systems with $SU(N)$ symmetry. In these
situations, the generators of $\mathfrak{su}(N)$ Lie algebra [92, 93] can be
applied. Another pathway to establish the theory of effective spectral density
function could be achieved by using the generalized Schwinger’s theory of
angular momentum, 333J. Schwinger, in Quantum Theory of Angular Momentum,
edited by L. C. Biedenharn and H. Van Dam (Academic, New York, 1965); see also
J. J. Sakurai, Modern Quantum Mechanics (Addison-Wesley, New York, 1994), p.
217. which remains to be further explored.

##  Appendix B The numerical accuracy of $t$-PFD strategies

In this section, we will provide examples of different $t$-PFD strategies to
illustrate its numerical accuracy on TCF fitting.

Fig. 5 presents the $t$-PFD results for the real and imaginary parts of bare-
bath TCFs. We take the most challenging model with bath parameters
$\omega_{c}/\Delta=1,\ \alpha=10$, and the bare-bath TCF plateau time is taken
as $40\Delta$ with resolution $dt=0.01\Delta$ uniformly. [62] Fig. 5a and b
present the accuracy of fitting for the zero-temperature model using different
number of terms. As is seen that using 4 or 5 terms will be accurate enough to
fit the real and imaginary parts of the TCF, respectively. Fig. 5c and d
present the fitting results of finite-temperature models using $4+4$ scheme.
As is observed straightforwardly that the real part of TCF is temperature-
independent, only the imaginary part varies with temperature, in accordance
with Eq. (23). The finite-temperature models are expected generally more
difficult to be accurately fitted than the zero-temperature model. Here all
curves are accurately fitted using $4+4$ scheme by eye inspection.

In summary, to reach better accuracy, one will need to use more terms, but the
expense is that the computational cost grow drastically. For practical use,
one will need to explore the proper $t$-PFD strategy to balance accuracy and
computational cost with regarding to the specific bath types and parameters.

## References

  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * Weiss [1993] U. Weiss,  _Quantum Dissipation Systems_ (World Scientific, Singapore, 1993).
  * Nitzan [2006] A. Nitzan,  _Chemical Dynamics in Condensed Phases_ (Oxford, New York, 2006).
  * Hänggi  _et al._ [1990] P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
  * Kofman and Kurizki [2004] A. G. Kofman and G. Kurizki, Phys. Rev. Lett. 93, 130406 (2004).
  * Hammes-Schiffer and Stuchebrukhov [2010] S. Hammes-Schiffer and A. A. Stuchebrukhov, Chemical Reviews 110, 6939 (2010).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Annals of Physics 149, 374 (1983).
  * Caldeira  _et al._ [1993] A. O. Caldeira, A. H. Castro Neto, and T. Oliveira de Carvalho, Phys. Rev. B 48, 13974 (1993).
  * Yan [2014] Y. J. Yan, Journal of Chemical Physics 140 (2014), 10.1063/1.4863379.
  * Prokof’ev and Stamp [1998] N. V. Prokof’ev and P. C. E. Stamp, Phys. Rev. Lett. 80, 5794 (1998).
  * Wernsdorfer  _et al._ [1999a] W. Wernsdorfer, T. Ohm, C. Sangregorio, R. Sessoli, D. Mailly, and C. Paulsen, Phys. Rev. Lett. 82, 3903 (1999a).
  * Wernsdorfer  _et al._ [1999b] W. Wernsdorfer, R. Sessoli, and D. Gatteschi, Europhysics Letters (EPL) 47, 254 (1999b).
  * Wernsdorfer  _et al._ [2000] W. Wernsdorfer, A. Caneschi, R. Sessoli, D. Gatteschi, A. Cornia, V. Villar, and C. Paulsen, Phys. Rev. Lett. 84, 2965 (2000).
  * Sinitsyn and Prokof’ev [2003] N. A. Sinitsyn and N. Prokof’ev, Phys. Rev. B 67, 134403 (2003).
  * Steane [1998] A. Steane, Reports on Progress in Physics 61, 117 (1998).
  * Dutt  _et al._ [2005] M. V. G. Dutt, J. Cheng, B. Li, X. Xu, X. Li, P. R. Berman, D. G. Steel, A. S. Bracker, D. Gammon, S. E. Economou, R.-B. Liu, and L. J. Sham, Phys. Rev. Lett. 94, 227403 (2005).
  * Rønnow  _et al._ [2005] H. M. Rønnow, R. Parthasarathy, J. Jensen, G. Aeppli, T. F. Rosenbaum, and D. F. McMorrow, Science 308, 389 (2005).
  * Du  _et al._ [2009] J. Du, X. Rong, N. Zhao, Y. Wang, J. Yang, and R. B. Liu, Nature 461, 1265 (2009).
  * Zhao  _et al._ [2012] N. Zhao, S.-W. Ho, and R.-B. Liu, Phys. Rev. B 85, 115303 (2012).
  * Zurek [1991] W. H. Zurek, Phys Today 37 (1991).
  * Anglin  _et al._ [1997] J. R. Anglin, J. P. Paz, and W. H. Zurek, Phys. Rev. A 55, 4041 (1997).
  * Álvarez  _et al._ [2015] G. A. Álvarez, D. Suter, and R. Kaiser, Science 349, 846 (2015).
  * Gómez-León and Stamp [2017] A. Gómez-León and P. C. E. Stamp, Phys. Rev. B 95, 054402 (2017).
  * Fay  _et al._ [2019] T. P. Fay, L. P. Lindoy, and D. E. Manolopoulos, The Journal of Chemical Physics 151, 154117 (2019).
  * Fay  _et al._ [2020] T. P. Fay, L. P. Lindoy, D. E. Manolopoulos, and P. J. Hore, Faraday Discuss. 221, 77 (2020).
  * Harvey and Wasielewski [2021] S. M. Harvey and M. R. Wasielewski, Journal of the American Chemical Society 143, 15508 (2021).
  * Bayliss  _et al._ [2020] S. L. Bayliss, D. W. Laorenza, P. J. Mintun, B. D. Kovos, D. E. Freedman, and D. D. Awschalom, Science 370, 1309 (2020).
  * Suárez and Silbey [1991] A. Suárez and R. Silbey, The Journal of Chemical Physics 95, 9115 (1991).
  * Makri [1999a] N. Makri, The Journal of Chemical Physics 111, 6164 (1999a).
  * Yan [2016] Y.-A. Yan, The Journal of Chemical Physics 145, 204111 (2016).
  * Makri [1999b] N. Makri, The Journal of Physical Chemistry B 103, 2823 (1999b).
  * Makri [2021] N. Makri, Phys. Chem. Chem. Phys. 23, 12537 (2021).
  * Wang and Shao [2012] H. Wang and J. Shao, The Journal of Chemical Physics 137, 22A504 (2012).
  * Gelman  _et al._ [2004] D. Gelman, C. P. Koch, and R. Kosloff, The Journal of Chemical Physics 121, 661 (2004).
  * Habecker  _et al._ [2019] F. Habecker, R. Röhse, and T. Klüner, The Journal of Chemical Physics 151, 134113 (2019).
  * Wu and Zhao [2013] N. Wu and Y. Zhao, The Journal of Chemical Physics 139, 054118 (2013).
  * Barnes  _et al._ [2011] E. Barnes, L. Cywiński, and S. Das Sarma, Phys. Rev. B 84, 155315 (2011).
  * Hsieh and Cao [2018a] C.-Y. Hsieh and J. Cao, The Journal of Chemical Physics 148, 014103 (2018a).
  * Hsieh and Cao [2018b] C.-Y. Hsieh and J. Cao, The Journal of Chemical Physics 148, 014104 (2018b).
  * Lü and Zheng [2009] Z. Lü and H. Zheng, The Journal of Chemical Physics 131, 134503 (2009).
  * Golosov  _et al._ [1999] A. A. Golosov, S. I. Tsonchev, P. Pechukas, and R. A. Friesner, The Journal of Chemical Physics 111, 9918 (1999).
  * Witzel and Das Sarma [2007] W. M. Witzel and S. Das Sarma, Phys. Rev. B 76, 241303 (2007).
  * Yao  _et al._ [2006] W. Yao, R.-B. Liu, and L. J. Sham, Phys. Rev. B 74, 195301 (2006).
  * Yang and Liu [2008] W. Yang and R.-B. Liu, Phys. Rev. B 78, 085315 (2008).
  * Yang  _et al._ [2020] Z.-S. Yang, Y.-X. Wang, M.-J. Tao, W. Yang, M. Zhang, Q. Ai, and F.-G. Deng, Annals of Physics 413, 168063 (2020).
  * Saikin  _et al._ [2007] S. K. Saikin, W. Yao, and L. J. Sham, Phys. Rev. B 75, 125314 (2007).
  * Gräßer  _et al._ [2021] T. Gräßer, P. Bleicker, D.-B. Hering, M. Yarmohammadi, and G. S. Uhrig, Phys. Rev. Research 3, 043168 (2021).
  * Nakamura and Tanimura [2021] K. Nakamura and Y. Tanimura, The Journal of Chemical Physics 155, 244109 (2021).
  * Kleinert [2009] H. Kleinert,  _Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets,5th ed._ (World Scientific, Singapore, 2009).
  * Lindoy and Manolopoulos [2018] L. P. Lindoy and D. E. Manolopoulos, Phys. Rev. Lett. 120, 220604 (2018).
  * Fay  _et al._ [2018] T. P. Fay, L. P. Lindoy, and D. E. Manolopoulos, The Journal of Chemical Physics 149, 064107 (2018).
  * Lindoy  _et al._ [2020] L. P. Lindoy, T. P. Fay, and D. E. Manolopoulos, The Journal of Chemical Physics 152, 164107 (2020).
  * Fay  _et al._ [2021] T. P. Fay, L. P. Lindoy, and D. E. Manolopoulos, The Journal of Chemical Physics 154, 084121 (2021).
  * Makarov and Makri [1994] D. E. Makarov and N. Makri, Chemical Physics Letters 221, 482 (1994).
  * Makri [1995] N. Makri, Journal of Mathematical Physics 36, 2430 (1995).
  * Segal  _et al._ [2010] D. Segal, A. J. Millis, and D. R. Reichman, Phys. Rev. B 82, 205323 (2010).
  * Meyer  _et al._ [2009] H. D. Meyer, F. Gatti, and G. A. Worth,  _Multidimensional Quantum Dynamics: MCTDH Theory and Applications_ (Wiley-VCH Verlag GmbH & Co. KGaA, 2009).
  * Tanimura [1990] Y. Tanimura, Phys. Rev. A 41, 6676 (1990).
  * Tanimura [2006] Y. Tanimura, Journal of the Physical Society of Japan 75, 082001 (2006).
  * Xu  _et al._ [2005] R.-X. Xu, P. Cui, X.-Q. Li, Y. Mo, and Y. Yan, The Journal of Chemical Physics 122, 041103 (2005).
  * Yan  _et al._ [2016] Y. Yan, J. Jin, R.-X. Xu, and X. Zheng, Frontiers of Physics 11, 110306 (2016).
  * Chen  _et al._ [2022] Z.-H. Chen, Y. Wang, X. Zheng, R.-X. Xu, and Y. Yan, The Journal of Chemical Physics 156, 221102 (2022).
  * Ising [1925] E. Ising, Zeitschrift fur Physik A Hadrons and Nuclei 31, 253 (1925).
  * Dirac [1929] P. A. M. Dirac, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 123, 714 (1929).
  * Note [1] There is a straightforward example with respect to the pure dephasing case, in which $\Delta=0$ such that the dephasing dynamics has analytical solution, as discussed by Rao and Kurizki [95], Hsieh and Cao [39].
  * Note [2] Note that historically, Schotte[96] first pointed out that at low temperature and long time limit ($t\gg\omega_{c}^{-1}$), the dynamics given by the Kondo Hamiltonian in a bosonic picture is the same as that of the corresponding fermionic operators. This pattern also emerges in the theoretical validation of the Caldeira-Leggett model [7, 8].
  * Burden and Faires [2011] R. L. Burden and J. D. Faires,  _Numerical Analysis, $9^{th}$ Ed._ (Brooks/Cole, Cengage Learning, 2011).
  * Shi  _et al._ [2009] Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan, The Journal of Chemical Physics 130, 084105 (2009).
  * Ozaki [2007] T. Ozaki, Phys. Rev. B 75, 035123 (2007).
  * Hu  _et al._ [2010] J. Hu, R.-X. Xu, and Y. Yan, The Journal of Chemical Physics 133, 101106 (2010).
  * Hu  _et al._ [2011] J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y. Yan, The Journal of Chemical Physics 134, 244106 (2011).
  * Liu  _et al._ [2014] H. Liu, L. Zhu, S. Bai, and Q. Shi, The Journal of Chemical Physics 140, 134106 (2014).
  * Duan  _et al._ [2017a] C. Duan, Q. Wang, Z. Tang, and J. Wu, The Journal of Chemical Physics 147, 164112 (2017a).
  * Wang  _et al._ [2019] Q. Wang, Z. Gong, C. Duan, Z. Tang, and J. Wu, The Journal of Chemical Physics 150, 084114 (2019).
  * Lambert  _et al._ [2019] N. Lambert, S. Ahmed, M. Cirio, and F. Nori, Nature Communications 10, 3721 (2019).
  * Anders  _et al._ [2007] F. B. Anders, R. Bulla, and M. Vojta, Phys. Rev. Lett. 98, 210402 (2007).
  * Wang and Thoss [2008] H. Wang and M. Thoss, New Journal of Physics 10, 115005 (2008).
  * Wang and Thoss [2010] H. Wang and M. Thoss, Chemical Physics 370, 78 (2010), dynamics of molecular systems: From quantum to classical.
  * Duan  _et al._ [2017b] C. Duan, Z. Tang, J. Cao, and J. Wu, Phys. Rev. B 95, 214308 (2017b).
  * Shao and Hänggi [1998] J. Shao and P. Hänggi, Phys. Rev. Lett. 81, 5710 (1998).
  * Forsythe and Makri [1999] K. M. Forsythe and N. Makri, Phys. Rev. B 60, 972 (1999).
  * Wilner  _et al._ [2015] E. Y. Wilner, H. Wang, M. Thoss, and E. Rabani, Phys. Rev. B 92, 195143 (2015).
  * Wang and Thoss [2006] H. Wang and M. Thoss, The Journal of Chemical Physics 124, 034114 (2006).
  * Jordan and Wigner [1929] P. Jordan and E. Wigner, Zeitschrift für Physik 47, 631 (1929).
  * Girardeau [1960] M. Girardeau, Journal of Mathematical Physics 1, 516 (1960).
  * Hall [2015] B. C. Hall,  _GTM222: Lie Groups, Lie Algebras, and Representations, An Elementary Introduction, $2^{nd}$ Ed._ (Springer, Switzerland, 2015).
  * Mattis [1965] D. C. Mattis,  _The Theory of Magnetism_ (Springer, 1965).
  * Wang  _et al._ [1966] Y. Wang, S. Shtrikman, and H. Callen, Journal of Applied Physics 37, 1451 (1966).
  * Jin  _et al._ [2007] J. Jin, S. Welack, J. Luo, X.-Q. Li, P. Cui, R.-X. Xu, and Y. Yan, The Journal of Chemical Physics 126, 134113 (2007).
  * Ghosh  _et al._ [2012] A. Ghosh, S. S. Sinha, and D. S. Ray, Phys. Rev. E 86, 011122 (2012).
  * Riera-Campeny  _et al._ [2021] A. Riera-Campeny, A. Sanpera, and P. Strasberg, PRX Quantum 2, 010340 (2021).
  * Runeson and Richardson [2020] J. E. Runeson and J. O. Richardson, The Journal of Chemical Physics 152, 084110 (2020).
  * Bossion  _et al._ [2022] D. Bossion, W. Ying, S. N. Chowdhury, and P. Huo, The Journal of Chemical Physics 157, 084105 (2022).
  * Note [3] J. Schwinger, in Quantum Theory of Angular Momentum, edited by L. C. Biedenharn and H. Van Dam (Academic, New York, 1965); see also J. J. Sakurai, Modern Quantum Mechanics (Addison-Wesley, New York, 1994), p. 217.
  * Bhaktavatsala Rao and Kurizki [2011] D. D. Bhaktavatsala Rao and G. Kurizki, Phys. Rev. A 83, 032105 (2011).
  * Schotte [1970] K. D. Schotte, Zeitschrift für Physik A Hadrons and nuclei 235, 155 (1970).