Simulating Open Quantum Systems with Hamiltonian Ensembles and the Nonclassicality of the Dynamics

TL;DR

This paper introduces a method to determine the nonclassical nature of open quantum system dynamics by analyzing whether they can be simulated with Hamiltonian ensembles, linking simulation impossibility to quantum correlations.

Contribution

It proposes a novel approach to assess system-environment correlations through Hamiltonian-ensemble simulation, defining nonclassicality based on simulation feasibility.

Findings

Nonunital dynamics are classified as nonclassical.

Examples include unital classical and unital nonclassical dynamics.

Simulation impossibility indicates nonclassical correlations.

Abstract

The incoherent dynamical properties of open quantum systems are generically attributed to an ongoing correlation between the system and its environment. Here, we propose a novel way to assess the nature of these system-environment correlations by examining the system dynamics alone. Our approach is based on the possibility or impossibility to simulate open-system dynamics with Hamiltonian ensembles. As we show, such (im)possibility to simulate is closely linked to the system-environment correlations. We thus define the nonclassicality of open-system dynamics in terms of the nonexistence of a Hamiltonian-ensemble simulation. This classifies any nonunital open-system dynamics as nonclassical. We give examples for open-system dynamics that are unital and classical, as well as unital and nonclassical.

Equations

$$\overline{\rho }(t)=\sum _j{p}_j\exp [-i{\hat{H}}_{j}t/\hslash ]{\rho }_0\exp [i{\hat{H}}_{j}t/\hslash ].$$

$$\frac{\partial }{\partial t}\overline{\rho }(t)=-\frac{i}{\hslash }[\epsilon (t){\hat{\sigma }}_z,\overline{\rho }(t)]+\gamma (t)[{\hat{\sigma }}_z\overline{\rho }(t){\hat{\sigma }}_z-\overline{\rho }(t)],$$

$$\varphi (t)={\int }_{-\infty }^{\infty }p(\omega )e^{i\omega t}d\omega .$$

$${\rho }_{\mathrm{T}}(t)=\sum _{k,j}{p}_j{\hat{E}}_{k,j}{\rho }_{\mathrm{S},0}{\hat{E}}_{k,j}^{†}\otimes |k⟩⟨k|,$$

$$\hat{U}(t)=\sum _j{\hat{E}}_j(t)\otimes \hat{\mathcal{U}}(t)|j⟩⟨j|.$$

$${\rho }_{\mathrm{T}}(t)=\sum _j{p}_j{e}^{-i{\hat{H}}_{j}t/\hslash }{\rho }_{\mathrm{S},0}{e}^{i{\hat{H}}_{j}t/\hslash }\otimes |j⟩⟨j|.$$

$${\hat{H}}_{\mathrm{S}}=\frac{\hslash {\omega }_0}{2}{\hat{\sigma }}_z$$

$${\hat{H}}_{\mathrm{I}}={\hat{\sigma }}_z\otimes$$

$$\varphi (t)=\exp \left[i{\omega }_{0}t-\Phi (t)\right].$$

$$\wp (\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\exp [i{\omega }_{0}t-\Phi (t)]e^{-i\omega t}dt.$$

$${\varphi }^{(\mathrm{X})}(t)=\exp \left[-i{\vartheta }_{\phi }(t)-\Phi (t)\right],$$

$${\vartheta }_{\phi }(t)$$

$${\rho }_{\mathrm{T}}(t)=\sum _{k,j}{p}_j{\hat{E}}_{k,j}{\rho }_{\mathrm{S},0}{\hat{E}}_{k,j}^{†}\otimes |k⟩⟨k|,$$

$$\sum _j{p}_j{\hat{E}}_{k,j}{\rho }_{\mathrm{S},0}{\hat{E}}_{k^{\prime },j}^{†}=0$$

$$\sum _k{\hat{E}}_{k,j}^{†}{\hat{E}}_{k,j}={\hat{E}}_{k_j,j}^{†}{\hat{E}}_{k_j,j}=\hat{I},\forall j.$$

$$\hat{U}(t)=\sum _j{\hat{E}}_j(t)\otimes \hat{\mathcal{U}}(t)|j⟩⟨j|.$$

$$\hat{U}(t+\delta t)=\hat{U}(t)\hat{U}(\delta t)$$

$$\hat{U}(t+\delta t)=\sum _j$$

$$\hat{U}(t)\hat{U}(\delta t)=\sum _{j^{\prime },j}$$

$$\frac{\partial \hat{\mathcal{U}}(t)}{\partial t}=\hat{\mathcal{U}}(t)\left[-\frac{i}{\hslash }\frac{\partial \hat{L}(0)}{\partial t}\right].$$

$${\rho }_{\mathrm{T}}(t)=\sum _j{p}_j{\hat{U}}_j{\rho }_{\mathrm{S},0}{\hat{U}}_j^{†}\otimes |j⟩⟨j|,$$

$$\sum _{j,k}f(t_j-t_k)z_j{z}_k^{\ast }\ge 0$$

$$\sum _{j,k}\varphi (t_j-t_k)z_j{z}_k^{\ast }=$$

$$\sum _{j,k}\exp \left[-\Phi (t_j-t_k)\right]\left(e^{i{\omega }_0{t}_j}{z}_j\right){\left(e^{i{\omega }_0{t}_k}{z}_k\right)}^{\ast }.$$

$$\mathcal{M}^{(2)}=\left[\begin{array}[]{cc}1&\exp\left[-\Phi(t_{2}-t_{1})\right]\\ \exp\left[-\Phi(t_{1}-t_{2})\right]&1\end{array}\right].$$

$$\mathcal{M}^{(3)}=\left[\begin{array}[]{ccc}1&e^{-\Phi_{2,1}}&e^{-\Phi_{3,1}}\\ e^{-\Phi_{1,2}}&1&e^{-\Phi_{3,2}}\\ e^{-\Phi_{1,3}}&e^{-\Phi_{2,3}}&1\end{array}\right],$$

$$\det \left({\mathcal{M}}^{(3)}\right)$$

$$\begin{array}[]{ccc}|z_{1}|^{2}&e^{-\Phi_{2,1}}z_{2}z_{1}^{\ast}&e^{-\Phi_{3,1}}z_{3}z_{1}^{\ast}\\ e^{-\Phi_{1,2}}z_{1}z_{2}^{\ast}&|z_{2}|^{2}&e^{-\Phi_{3,2}}z_{3}z_{2}^{\ast}\\ e^{-\Phi_{1,3}}z_{1}z_{3}^{\ast}&e^{-\Phi_{2,3}}z_{2}z_{3}^{\ast}&|z_{3}|^{2}\end{array}.$$

$$\begin{array}[]{cccc}|z_{1}|^{2}&e^{-\Phi_{2,1}}z_{2}z_{1}^{\ast}&e^{-\Phi_{3,1}}z_{3}z_{1}^{\ast}&e^{-\Phi_{4,1}}z_{4}z_{1}^{\ast}\\ e^{-\Phi_{1,2}}z_{1}z_{2}^{\ast}&|z_{2}|^{2}&e^{-\Phi_{3,2}}z_{3}z_{2}^{\ast}&e^{-\Phi_{4,2}}z_{4}z_{2}^{\ast}\\ e^{-\Phi_{1,3}}z_{1}z_{3}^{\ast}&e^{-\Phi_{2,3}}z_{2}z_{3}^{\ast}&|z_{3}|^{2}&e^{-\Phi_{4,3}}z_{4}z_{3}^{\ast}\\ e^{-\Phi_{1,4}}z_{1}z_{4}^{\ast}&e^{-\Phi_{2,4}}z_{2}z_{4}^{\ast}&e^{-\Phi_{3,4}}z_{3}z_{4}^{\ast}&|z_{4}|^{2}\end{array},$$

$$\begin{array}[]{cc}|z_{1}+z_{2}+z_{3}|^{2}&e^{-\Phi_{4,1}}z_{4}(z_{1}+z_{2}+z_{3})^{\ast}\\ e^{-\Phi_{1,4}}(z_{1}+z_{2}+z_{3})z_{4}^{\ast}&|z_{4}|^{2}\end{array}.$$

$${\hat{H}}_{\mathrm{S}}$$

Full text

Introduction.—When a quantum system interacts with its environment, its dynamical behavior will, in general, deviate from the dynamics of a strictly isolated one Leggett et al. (1987); Breuer and Petruccione (2002); Weiss (2012); Zhang et al. (2012); Xiong et al. (2015); Chen et al. (2015b). As a result of an ongoing bipartite correlation arising from the system-environment interaction, the system dynamics may display incoherent characteristics, such as dephasing or damping processes. Formally, such processes are captured by quantum master equations, replacing the von Neumann equation for isolated systems.

However, incoherent dynamics can also arise as a consequence of a purely classical averaging procedure over distinct autonomous evolutions. For example, the double slit experiment can, when exposed to a disordered potential and after averaging, encounter similar decoherence as if which-slit information had leaked into an environment Gneiting et al. (2016). In this sense, disordered quantum systems described by Hamiltonian ensembles can behave in an analogous manner as open quantum systems, even if individual realizations are strictly isolated (Fig. 1) Gneiting et al. (2016); Kropf et al. (2016); Gneiting and Nori (2017, 2017b).

Here, we exploit this dynamical correspondence to assess the nature of the system-environment correlations in terms of the system properties alone. As we show, the impossibility to simulate is necessarily linked to nonclassical system-environment correlations. On the other hand, if such a simulation is possible, then there always exists a system-environment model which reproduces the system dynamics by relying only on classical correlations. This leads us to defining the nonclassicality of open-system dynamics in terms of the nonexistence of a Hamiltonian-ensemble simulation.

Alternative definitions for the nonclassicality of system dynamics have been proposed Rahimi-Keshari et al. (2013); Sabapathy (2016). In these definitions, the dynamics is considered classical if the state preserves classicality during the temporal evolution. Typically, the classicality of states in these approaches is formulated in terms of the Wigner function or the Glauber-Sudarshan $P$ representation Wigner (1932); Glauber (1963); Sudarshan (1963); Miranowicz et al. (2010); Bartkowiak et al. (2011); Miranowicz et al. (2015). While these definitions also rely on system properties alone, their applicability is limited to systems amenable for such a phase space description, excluding other cases of interest. We, instead, propose to discuss the nonclassicality of open-system dynamics separately, based on the system-environment correlations and independent of the nature of the system.

As an immediate consequence of our definition, all nonunital dynamics, e.g., dissipative processes, are nonclassical. The spin-boson model, in contrast, which displays unital dynamics and—on the level of the model—quantum correlations, can, as we find, be simulated and thus exhibits classical open-system dynamics according to our definition, since these correlations cannot be certified by considering the system dynamics alone. In the case of an extended spin-boson model, however, where the environment is complemented by a second qubit and the system dynamics remains unital, we prove the nonexistence of a simulating Hamiltonian ensemble for certain spectral densities; i.e., the system’s evolution is in these cases witnessed to be manifestly nonclassical.

Dynamics of Hamiltonian ensembles.—An isolated quantum system is described by a Hamiltonian ensemble (HE) $\{(p_{j},\widehat{H}_{j})\}$, if the autonomous Hamiltonian $\widehat{H}_{j}$ of the system is drawn from a probability distribution $p_{j}$ [see Fig. 1(a)]. Such HEs are applicable to describing disordered quantum systems. Here, we relate HEs to open quantum systems.

The dynamics of the ensemble averaged state $\overline{\rho}(t)$ exhibits features distinct from the dynamics of any single realization. The latter is governed by the unitary evolution $\rho_{j}(t)=\widehat{U}_{j}\rho_{0}\widehat{U}_{j}^{\dagger}$, with the initial state $\rho_{0}$ and $\widehat{U}_{j}=\exp[-i\widehat{H}_{j}t/\hbar]$, whereas the dynamics of the averaged state $\overline{\rho}(t)$ is given by the unital (i.e., identity invariant) map

[TABLE]

Note that an evolution equation for $\overline{\rho}(t)$ cannot be reduced to some effective Hamiltonian alone but must, in general, take the form of a quantum master equation Gneiting et al. (2016); Kropf et al. (2016) (see Xu et al. (2013) for an experimental implementation of this).

A seminal and instructive example considers a single qubit subject to spectral disorder; i.e., the Hamiltonians in the ensemble differ only in their eigenvalues, while they share a common basis of eigenstates Kropf et al. (2016). The HE may be given by $\{(p(\omega),\hbar\omega\hat{\sigma}_{z}/2)\}$, with the probability distribution $p(\omega)$ kept general. The resulting master equation reads

[TABLE]

where the effective energy $\varepsilon(t)=\hbar\mathrm{Im}[\partial_{t}\ln\phi(t)]/2$ and the decoherence rate $\gamma(t)=-\mathrm{Re}[\partial_{t}\ln\phi(t)]/2$ follow from the dephasing factor

[TABLE]

Depending on the underlying probability distribution $p(\omega)$, the master equation (2) can range from time-constant dephasing to a strongly oscillating incoherent behavior, the latter even giving rise to purity revivals Kropf et al. (2016).

It is worthwhile to recall that the occurrence of incoherent dynamics in the case of HEs is a consequence of the averaging procedure. Nevertheless, it is reminiscent of open quantum systems, where, in contrast, an ongoing correlation between the system and environment gives rise to the incoherent dynamics. This Letter explores the possibility to simulate open quantum systems with HEs, and vice versa, and the implications on the system-environment correlations.

Simulating open quantum systems with Hamiltonian ensembles.—We now show that nonclassical system-environment correlations are necessarily linked to the impossibility to simulate the open-system dynamics with a HE. To this end, we show—conversely—that, if system and environment are persistently classically correlated, then the reduced system state is described by a HE.

A system-environment arrangement is characterized by an autonomous total Hamiltonian $\widehat{H}_{\mathrm{T}}=\widehat{H}_{\mathrm{S}}+\widehat{H}_{\mathrm{E}}+\widehat{H}_{\mathrm{I}}$, with the system $\widehat{H}_{\mathrm{S}}$, the environment $\widehat{H}_{\mathrm{E}}$, and the interaction Hamiltonian $\widehat{H}_{\mathrm{I}}$ [see Fig. 1(b)]. The total system evolves unitarily as $\rho_{\mathrm{T}}(t)=\widehat{U}\rho_{\mathrm{T},0}\widehat{U}^{\dagger}$, with $\widehat{U}=\exp[-i\widehat{H}_{\mathrm{T}}t/\hbar]$. We say that an open system is described by a HE if the reduced system state $\rho_{\mathrm{S}}(t)=\mathrm{Tr}_{\mathrm{E}}[\rho_{\mathrm{T}}(t)]$ allows a decomposition of the form (1), where the probabilities $p_{j}$ and the Hamiltonians $\widehat{H}_{j}$ of the ensemble are determined by $\widehat{H}_{\mathrm{T}}$ and the initial state $\rho_{\mathrm{T},0}$.

Instead of further specifying the total Hamiltonian $\widehat{H}_{\mathrm{T}}$, we now assume that the total state $\rho_{\mathrm{T}}(t)$ remains at all times classically correlated, displaying neither quantum discord Ollivier and Zurek (2001); Dakić et al. (2010) nor entanglement. Under this condition, we argue that, for every initial state of the form $\rho_{\mathrm{T},0}=\rho_{\mathrm{S},0}\otimes\sum_{j}p_{j}{|j\rangle}{\langle j|}$ (with $\{p_{j}\}$ a time-independent probability distribution and $\{{|j\rangle}\}$ a basis of the environment), the reduced system state can be described by a HE.

The detailed proof is presented in the Supplemental Material sup . Here, we outline the central steps. First, as a direct consequence of the classical correlations, there exists an environmental basis $\{{|k\rangle}\}$ (in general different from $\{{|j\rangle}\}$), such that

[TABLE]

where the operators $\widehat{E}_{k,j}={\langle k|}\widehat{U}(t){|j\rangle}$ act on the system and satisfy $\sum_{k}\widehat{E}_{k,j}^{\dagger}\widehat{E}_{k,j}=\widehat{I}$ for each $j$.

To demonstrate that the $\widehat{E}_{k,j}$ are unitary, we again use the zero-discord assumption, which implies that each environmental off-diagonal term vanishes, i.e., $\widehat{E}_{k,j}\rho_{\mathrm{S},0}\widehat{E}_{k^{\prime},j}^{\dagger}=0$ for $k\neq k^{\prime}$. This, in turn, implies that there exists a bijection between $\{{|j\rangle}\}$ and $\{{|k\rangle}\}$ such that $\widehat{E}_{k,j}$ is nonzero only when its two indices match the bijection, i.e., $\widehat{E}_{k,j}=\widehat{E}_{k_{j^{\prime}},j}\delta_{j,j^{\prime}}$. Unitarity of the $\widehat{E}_{k,j}$ then follows directly.

Finally, we address the time dependence of the $\widehat{E}_{k,j}$. Expressing the bijection as a unitary operator $\widehat{\mathcal{U}}(t)$ and safely neglecting the index $k$, we can recast $\widehat{U}(t)$ in a separable form:

[TABLE]

The group properties of $\{\widehat{U}(t)|t\in\mathbb{R}\}$ are thus inherited by the operators $\widehat{E}_{j}(t)$ and $\widehat{\mathcal{U}}(t)$; i.e., due to the time independence of the total Hamiltonian, we can write $\widehat{E}_{j}(t)=\exp[-i\widehat{H}_{j}t/\hbar]$ and $\widehat{\mathcal{U}}(t)=\sum_{j}\exp[-i(\theta_{j}t/\hbar)]{|j\rangle}{\langle j|}$, with $\widehat{H}_{j}$ time-independent Hermitian operators and $\theta_{j}$ real-valued constants. Consequently, Eq. (1) corresponds to a time-independent HE $\{(p_{j},\widehat{H}_{j})\}$ when tracing over the environment.

Simulating Hamiltonian ensembles with open quantum systems.—The impossibility to simulate an open system with a HE certifies the quantum nature of the system-environment correlations. Notably, this is achieved by considering system properties alone, i.e., without explicit reference to the environment. We now show that, on the other hand, the existence of a simulating HE always admits the possibility of classical system-environment correlations; i.e., the latter cannot be excluded by considering only the system.

We explicitly construct a system-environment arrangement which reproduces an arbitrary HE $\{(p_{j},\widehat{H}_{j})\}$ relying only on classical correlations. To this end, we write $\widehat{H}_{j}=\widehat{\overline{H}}+\widehat{V}_{j}$ (with the average $\widehat{\overline{H}}=\sum_{j}p_{j}\widehat{H}_{j}$) and choose the interaction to be of the form $\widehat{H}_{\rm I}=\sum_{j}\widehat{V}_{j}\otimes{|j\rangle}{\langle j|}$; i.e., it associates with each $\widehat{H}_{j}$ of the ensemble a distinct state ${|j\rangle}$ of an (arbitrary) basis of the environment. Note that the index $j$ is generic and may be continuous and/or a multi-index. The environment must then be chosen appropriately to accommodate the complexity of the HE. Moreover, we take the system Hamiltonian to be the average $\widehat{H}_{\rm S}=\widehat{\overline{H}}$ and the bath Hamiltonian to be diagonal in the same basis as $\widehat{H}_{\rm I}$, i.e., $[\widehat{H}_{\rm E},\widehat{H}_{\rm I}]=0$.

With a separable initial state, $\rho_{\rm T}(0)=\rho_{{\rm S},0}\otimes\rho_{\rm E}$, and $\rho_{\rm E}=\sum_{j}p_{j}{|j\rangle}{\langle j|}$ (i.e., $[\rho_{\rm E},\widehat{H}_{\rm E}]=0$, and the probabilities of the Hamiltonian ensemble are assigned to the environmental populations), the time-evolved total state reads $\rho_{\rm T}(t)=\widehat{U}_{\rm S+I}(\rho_{\rm S,0}\otimes\rho_{\rm E})\widehat{U}^{\dagger}_{\rm S+I}$, with $\widehat{U}_{\rm S+I}=\exp[-i(\widehat{H}_{\rm S}+\widehat{H}_{\rm I})t/\hbar]$. Rewriting $\widehat{U}_{\rm S+I}=\sum_{j}\widehat{U}_{j}\otimes{|j\rangle}{\langle j|}$, with $\widehat{U}_{j}=\exp[-i\widehat{H}_{j}t/\hbar]$, we obtain

[TABLE]

If we now trace over the environment, $\rho_{\rm S}(t)=\sum_{j}{\langle j|}\rho_{\rm T}(t){|j\rangle}$, we recover the desired decomposition (1) in terms of the HE $\{(p_{j},\widehat{H}_{j})\}$. Moreover, it is easy to see that the total state (6) is exclusively classically correlated, as desired.

As an instructive example, we consider a pair of qubits coupled to each other via a controlled-NOT gate, where a control (C) qubit determines the operation on a target (T) qubit. If the state of the C qubit resides in the classical mixture $\rho_{\mathrm{C}}=a{|1\rangle}{\langle 1|}+(1-a){|0\rangle}{\langle 0|}$ Chen et al. (2015a), the reduced dynamics of the T qubit will be described by the mixture of evolutions $\rho_{\mathrm{T}}(t)=a\widehat{U}_{x}\rho_{\mathrm{T},0}\widehat{U}^{\dagger}_{x}+(1-a)\rho_{\mathrm{T},0}$, with $\widehat{U}_{x}=\exp[-iJ\hat{\sigma}_{x}t/2\hbar]$ and $J$ the coupling strength. We thus recover the HE $\{(a,J\hat{\sigma}_{x}/2),(1-a,\widehat{I})\}$. In this example, the C qubit plays the role of environment, and the qubit pair is at most classically correlated.

Nonclassicality of the dynamics.—It appears natural to classify open-system dynamics according to their correlation with the environment; i.e., if the system and environment are merely classically correlated, the dynamics may be considered classical; if they are quantum correlated, one may call the dynamics nonclassical. In most cases, however, one does not have (full) access to the environment, rendering such an immediate definition problematic.

We now suggest to classify open-system dynamics by the (im)possibility to describe the system dynamics by a HE. On the one hand, this definition relies only on system properties, as desired from a practical point of view. On the other hand, as we have shown, it directly links to the system-environment correlations, as desirable from a conceptual perspective. Whenever such a simulation exists, it is impossible to exclude classical system-environment correlations by knowledge of the system dynamics alone, and we call the latter classical. If the simulation does not exist, quantum correlations must be involved; hence, the dynamics is nonclassical.

As a direct consequence of our definition, any nonunital dynamics is classified nonclassical—a simulating HE is manifestly excluded, certifying the presence of quantum correlations. This includes, e.g., dissipative processes such as the spontaneous decay of an atom. On the other hand, according to our operational definition, we may even call an open-system dynamics classical if the actual system-environment correlations are quantum. This is because our approach is deliberately ignorant of the actual environment and relies only on the possibility to explain the system dynamics with classical correlations. Next, we give an example for this.

Simulating the spin-boson model.—We now show that the system dynamics of the spin-boson model can be simulated by a HE, even though the actual model displays quantum correlations Bera et al. (2014); Iles-Smith et al. (2014). The spin-boson model

[TABLE]

has been extensively studied and is analytically solvable Breuer and Petruccione (2002). Tracing over the environment, the qubit system exhibits pure dephasing dynamics characterized by the dephasing factor

[TABLE]

In contrast to Eq. (3), which results from averaging over a HE, the dephasing factor (8) incorporates the information of the interaction and the environment into $\Phi(t)=4\int_{0}^{\infty}\omega^{-2}\mathcal{J}(\omega)\coth\left(\hbar\omega/2k_{\mathrm{B}}T\right)\left(1-\cos\omega t\right)d\omega$, where $\mathcal{J}(\omega)$ is the environmental spectral density. The above solution assumes that the initial state is a direct product, and that the environment is initially thermalized at temperature $T$.

We now construct a simulating HE. In view of Eq. (2), we deduce that individual member Hamiltonians in the ensemble must be of the form $\omega\hat{\sigma}_{\mathrm{z}}/2$, which leaves us with determining the corresponding probabilities. Given a probability distribution $p(\omega)$, the averaged dynamics can be determined by Eq. (3). Conversely, the underlying distribution function leading to a specific dephasing factor (8) is obtained via the inverse Fourier transform

[TABLE]

It is clear that the effect of $\omega_{0}$ is merely to shift $\wp(\omega)$.

To be a legitimate probability distribution function, the resulting $\wp(\omega)$ in Eq. (9) must be normalized [$\int_{-\infty}^{\infty}\wp(\omega)d\omega=1$], real [$\wp(\omega)\in\mathbb{R}$], and positive [$\wp(\omega)\geq 0$]. Normalization is easily seen, since $\phi(0)=1$ follows from the fact that the pure dephasing dynamics, characterized by Eq. (8), should be completely positive and trace preserving. We therefore have $\int_{-\infty}^{\infty}\wp(\omega)d\omega=(2\pi)^{-1}\int_{-\infty}^{\infty}\exp\left[i\omega_{0}t-\Phi(t)\right]2\pi\delta(t-0)dt=1$. Moreover, since one is generically interested in the dynamical properties only for $t\geq 0$, we can deliberately extend the time domain to the full real axis such that $\Phi(t)$ is even and $\phi(-t)=\phi(t)^{\ast}$. This guarantees that $\wp(\omega)$ is real: $\wp(\omega)=(\pi)^{-1}\int_{0}^{\infty}\exp\left[-\Phi(t)\right]\cos(\omega-\omega_{0})tdt\in\mathbb{R}$.

The positivity of $\wp(\omega)$ is less obvious, due to the sinusoidal factors of the integrand in Eq. (9). In the following, we invoke Bochner’s theory Bochner (1933) to prove the general positivity of $\wp(\omega)$. To this end, we first introduce the notion of positive definiteness. A function $f:\mathbb{R}\rightarrow\mathbb{C}$ is called positive definite if it satisfies $\sum_{j,k}f(t_{j}-t_{k})z_{j}z_{k}^{\ast}\geq 0$ for any finite number of pairs $\{(t_{j},z_{j})|t_{j}\in\mathbb{R},z_{j}\in\mathbb{C}\}$. Note that positive definiteness of a function is different from a positive function, since the latter may not necessarily be positive definite and vice versa. Rather, it corresponds to the positive semidefiniteness of a Hermitian matrix $\left[f(t_{j}-t_{k})\right]_{j,k\in\mathcal{S}}$, formed by the function values $f(t_{j}-t_{k})$ in accordance with a certain set of indices $\mathcal{S}$. As one can show, $\phi(t)$ in Eq. (8) is indeed positive definite. The proof is given in Ref. sup .

Bochner’s theorem states that a function $f$, defined on $\mathbb{R}$, is the Fourier transform of unique positive measure with density function $\wp$, if and only if $f$ is continuous and positive definite Loomis (1953); Rudin (1990). We can thus conclude that $\phi(t)$ in Eq. (8) is the Fourier transform of a certain valid probability distribution $\wp(\omega)$ [Eq. (9)], i.e., an analog to Eq. (3).

In summary, we have proven that there exists a unique HE, $\left\{\left(\wp(\omega),\omega\hat{\sigma}_{\mathrm{z}}/2\right)\right\}$, which simulates the system dynamics exactly, irrespective of the spectral density $\mathcal{J}(\omega)$ and the associated, possibly intricate system-environment entanglement. We thus call this dynamics classical.

Extended spin-boson model.—Unless the system dynamics is nonunital, proving the nonexistence of a simulating HE is, in general, a nontrivial task. We now accomplish this for an extended spin-boson model, at the same time deducing the presence of quantum correlations from system properties alone.

Our model consists of two qubits coupled to a common boson environment. The system and the interaction Hamiltonian are replaced by $\widehat{H}_{\mathrm{S}}=\sum_{j}\hbar\omega_{j}\hat{\sigma}_{z,j}/2$ ($j=1$, $2$) and $\widehat{H}_{\mathrm{I}}=\sum_{j,\vec{k}}\hat{\sigma}_{z,j}\otimes\hbar(g_{j,\vec{k}}\hat{b}_{\vec{k}}^{\dagger}+g_{j,\vec{k}}^{*}\hat{b}_{\vec{k}})$, respectively, while the environment Hamiltonian $\widehat{H}_{\mathrm{E}}$ is kept as in Eq. (Simulating Open Quantum Systems with Hamiltonian Ensembles and the Nonclassicality of the Dynamics). Note that the two qubits do not interact directly. The coupling constants $g_{j,\vec{k}}$ are, in general, complex numbers. In order to reveal the nonclassical effects caused by their relative phase, we assume, for simplicity, that they have the same amplitude, i.e., $g_{2,\vec{k}}=g_{1,\vec{k}}e^{i\varphi}$.

In the interaction picture, the total system evolves according to $\widehat{U}^{\mathrm{I}}(t)=\mathcal{T}\left\{\exp\left[-i\int_{0}^{t}\sum_{\vec{k}}\widehat{Z}_{\vec{k}}\hat{b}_{\vec{k}}^{\dagger}(\tau)+\widehat{Z}_{\vec{k}}^{\dagger}\hat{b}_{\vec{k}}(\tau)d\tau\right]\right\}$, with $\mathcal{T}$ the time-ordering operator, $\widehat{Z}_{\vec{k}}=\sum_{j}g_{j,\vec{k}}\hat{\sigma}_{z,j}$, and $\hat{b}_{\vec{k}}(t)=e^{-i\omega_{\vec{k}}t}\hat{b}_{\vec{k}}$, respectively. In contrast to the conventional spin-boson model, time ordering plays a nontrivial role here Reina et al. (2002). (For details, see Supplemental Material sup .)

In the following, we regard one qubit as the system and the other as part of the environment. The reduced dynamics of the system qubit is then pure dephasing with the dephasing factor [cf. Eq. (8)]

[TABLE]

where

[TABLE]

In the second line, we have manually inserted $\mathrm{sign}(t)$. This ensures that $\phi^{(\mathrm{X})}(-t)=\phi^{(\mathrm{X})*}(t)$ and $\wp^{(\mathrm{X})}(\omega)\in\mathbb{R}$. The presence of $\vartheta_{\varphi}(t)$, however, will, in general, result in the violation of positivity. Note that, similar to the conventional spin-boson model, individual member Hamiltonians in the HE must be of the form $\omega\hat{\sigma}_{\mathrm{z}}/2$, which allows us to follow the same line of argument.

To demonstrate this violation, we consider the Ohmic spectral density $\mathcal{J}_{\mathrm{o1}}(\omega)=\omega\exp(-\omega/\omega_{\mathrm{c}})$ in the zero-temperature limit and a degenerate system Hamiltonian, i.e., $\omega_{j}=0$. In Fig. 2(a), we depict the legitimate probability distribution $\wp_{\mathrm{o1}}(\omega)$ for the conventional spin-boson model at $\omega_{\mathrm{c}}=1$ (blue curve) and $\omega_{\mathrm{c}}=3$ (red curve), while in Fig. 2(b), we show $\wp_{\mathrm{o1}}^{(\mathrm{X})}(\omega)$ for our extended model with $\varphi=\pi/4$ (solid curves) and $\varphi=5\pi/4$ (dashed curves). The latter display a manifest violation of positivity. In Fig. 2(c), we show the landscape of negative contributions to $\wp_{\mathrm{o1}}^{(\mathrm{X})}(\omega)$ against $\omega$ and $\varphi$ for $\omega_{\mathrm{c}}=1$. The gray dashed lines highlight $\varphi=\pi/4$ and $5\pi/4$, chosen in Fig. 2(b).

Conclusions.—We propose a way to classify open-system dynamics according to their system-environment correlations, i.e., if the latter are classical or quantum. As we showed, this can be tested by knowledge of the system evolution alone, based on the (im)possibility to simulate the open-system dynamics with a Hamiltonian ensemble. According to our definition, any nonunital dynamics is nonclassical. Some unital system evolutions, however, such as in the spin-boson model, are classified as classical, even though the model displays quantum correlations. This highlights the operational nature of our definition.

With the extended spin-boson model, we provide an example for unital dynamics which is nonclassical according to our definition. Let us note that one may be able to simulate a larger class of unital dynamics with time-dependent Hamiltonian ensembles. It is, for example, known that, in the case of qubits, any unital dynamics can be simulated with an ensemble of time-dependent Hamiltonians, albeit only if also the probabilities are allowed to be time dependent Landau and Streater (1993); Audenaert and Scheel (2008). However, in the case of autonomous system-environment arrangements, which we consider here, such generalization appears unjustified. Finally, let us remark that demonstrating the nonexistence of a Hamiltonian-ensemble simulation is, in the case of unital evolutions, in general, a nontrivial task. An equivalent but simpler test appears desirable.

Acknowledgements.—This work is supported partially by the National Center for Theoretical Sciences and Ministry of Science and Technology, Taiwan, Grants No. MOST 103-2112-M-006-017-MY4 and No. MOST 105-2811-M-006-059, the MURI Center for Dynamic Magneto-Optics via AFOSR Award No. FA9550-14-1-0040, the Japan Society for the Promotion of Science (KAKENHI), the IMPACT program of JST, JSPS-RFBR Grant No. 17-52-50023, CREST Grant No. JPMJCR1676, RIKEN-AIST Challenge Research Fund, and the Sir John Templeton Foundation.

I TIME-INDEPENDENT HAMILTONIAN ENSEMBLE

In the following, we elaborate in detail the proof that classical bipartite correlations allow for a time-independent Hamiltonian ensemble decomposition of the reduced system dynamics. In this proof, we do not assume a specific form of the total Hamiltonian $\widehat{H}_{\mathrm{T}}$. To be precise, we now show that, if the total Hamiltonian $\widehat{H}_{\mathrm{T}}$ is time-independent, and if, for every initial state of the form $\rho_{\mathrm{T},0}=\rho_{\mathrm{S},0}\otimes\sum_{j}p_{j}{|j\rangle}{\langle j|}$ (with $\{p_{j}\}$ being any time-independent probability distribution), the time-evolved total state $\rho_{\mathrm{T}}(t)=\widehat{U}(t)\rho_{\mathrm{T},0}\widehat{U}^{\dagger}(t)$ is always classically correlated between system and environment, displaying neither quantum discord nor entanglement at any time, then the reduced system dynamics admits a time-independent Hamiltonian ensemble decomposition.

Proof.

Due to the zero-discord assumption, there exists an environmental basis $\{{|k\rangle}\}$ (in general time-dependent and different from $\{{|j\rangle}\}$), such that

[TABLE]

where $\widehat{E}_{k,j}={\langle k|}\widehat{U}(t){|j\rangle}$ are operators acting on the system Hilbert space satisfying $\sum_{k}\widehat{E}_{k,j}^{\dagger}\widehat{E}_{k,j}=\widehat{I}$ for each $j$. At this point, we are not yet clear about the time-dependence of $\widehat{E}_{k,j}$ and ${|k\rangle}$ nor the unitarity of $\widehat{E}_{k,j}$.

Crucially, the condition

[TABLE]

should hold for any $k\neq k^{\prime}$, due to the zero-discord assumption. Therefore each term $\widehat{E}_{k,j}\rho_{\mathrm{S},0}\widehat{E}_{k^{\prime},j}^{\dagger}$ in the above equation vanishes individually. The only possibility to reconcile Eqs. (1) and (2) is the existence of a specific bijection between $\{{|j\rangle}\}$ and $\{{|k\rangle}\}$, such that $\widehat{E}_{k,j}=\widehat{E}_{k_{j^{\prime}},j}\delta_{j,j^{\prime}}$ for each $j$, i.e., $\widehat{E}_{k,j}$ is non-zero only when its two indices match the bijection. Then the unitarity of $\widehat{E}_{k_{j},j}$ can then be confirmed according to

[TABLE]

The bijection between $\{|j\rangle\}$ and $\{|k\rangle\}$ can be expressed in terms of a unitary operator $\widehat{\mathcal{U}}(t)$, such that ${\langle k_{j^{\prime}}|}\widehat{\mathcal{U}}(t){|j\rangle}=\delta_{j,j^{\prime}}$. The unitary evolution operator can then be recast in a separable form,

[TABLE]

In the following discussion, we can, in order to keep the notation simple, safely neglect the index $k$.

Since $\{\widehat{U}(t)=\exp[-i\widehat{H}_{\mathrm{T}}t/\hbar]|t\in\mathbb{R}\}$ forms a group isomorphism on $\mathbb{R}$, we have the one-parameter group property

[TABLE]

for $t\in\mathbb{R}$ and infinitesimal $\delta t$. Due to the unitarity of $\widehat{\mathcal{U}}(t)$, it can be expressed in terms of an Hermitian generator $\widehat{L}(t)$ in the $\mathfrak{u}(\mathrm{dim}\mathcal{H}_{\mathrm{E}})$ Lie algebra on the environmental Hilbert space $\mathcal{H}_{\mathrm{E}}$ such that $\widehat{\mathcal{U}}(t)=\exp[-i\widehat{L}(t)/\hbar]$. Together with Eq. (4), the left hand side of Eq. (5) can be written as

[TABLE]

On the right hand side of Eq. (I), we expand $\widehat{\mathcal{U}}(t+\delta t)$ around $t$ to first order in $\delta t$. Notably, since we do not know the time-dependence and commutativity of $\widehat{L}(t)$ at this point, we can only achieve a formal expansion in Eq. (I).

Meanwhile, the right hand side of Eq. (5) reads

[TABLE]

We again expand $\widehat{\mathcal{U}}(\delta t)$ around $t=0$. However, unlike the formal expansion in Eq. (I), we now obtain an explicit expansion in Eq. (I), since $\widehat{\mathcal{U}}(0)=\widehat{I}$ commutes with any operator.

Comparing Eqs. (I) and (I), we conclude from their first terms that the group property $\widehat{E}_{j}(t+\delta t)=\widehat{E}_{j}(t)\widehat{E}_{j}(\delta t)$ holds and, combined with the unitarity inferred in Eq. (3), that time-independent Hermitian operators $\widehat{H}_{j}$ exist, such that $\widehat{E}_{j}(t)=\exp[-i\widehat{H}_{j}t/\hbar]$, as well.

To reconcile the second terms of Eqs. (I) and (I), $\partial\widehat{L}(0)/\partial t$ should be diagonalized in the basis $\{{|j\rangle}\}$, such that $\partial\widehat{L}(0)/\partial t=\sum_{j}(\partial\theta_{j}(0)/\partial t){|j\rangle}{\langle j|}$, with real parameters $\theta_{j}(t)$. Moreover, $\widehat{\mathcal{U}}(t)$ should satisfy

[TABLE]

To guarantee its validity, $\partial\widehat{L}(0)/\partial t$ should commute with $\widehat{L}(t)$, since the latter is the generator of $\widehat{\mathcal{U}}(t)$. Consequently, the time-dependence of each $\theta_{j}(t)$ can be of first order, such that $\widehat{\mathcal{U}}(t)=\sum_{j}\exp[-i(\theta_{j}t/\hbar)]{|j\rangle}{\langle j|}$, with real constants $\theta_{j}$. ∎

Consequently, the total state in Eq. (1) can be rewritten as

[TABLE]

with $\widehat{U}_{j}=\exp[-i\widehat{H}_{j}t/\hbar]$, which corresponds to a time-independent Hamiltonian ensemble $\{(p_{j},\widehat{H}_{j})\}$ when tracing over the environment.

Finally, let us remark that, while we restrict ourselves to a time-independent total Hamiltonian, some of our conclusions can be easily generalized to the time-dependent case. This is because Eqs. (1-4) are consequences of the zero-discord assumption alone, regardless of the time-dependence of the total Hamiltonian. Therefore, we can also achieve the ensemble form with time-varying member Hamiltonians for a time-dependent total Hamiltonian. However, as discussed in the main article, in the case of autonomous system-environment arrangements, i.e., in the absence of external control, such generalization appears unjustified.

Additionally, we note that, for the case of time-independent total Hamiltonians, the separable form (4) not only guarantees a persistently classically correlated total state, but also keeps the environmental basis intact without rotation, up to a phase angle $\theta_{j}t$.

II POSITIVE DEFINITENESS

Here, we present the proof of the positive definiteness of the dephasing factor $\phi(t)=\exp\left[i\omega_{0}t-\Phi(t)\right]$. For completeness, we recall the definition of positive definiteness.

Positive definiteness: A function $f$ defined on $\mathbb{R}$ is called positive definite if it satisfies

[TABLE]

for any finite number of pairs $\{(t_{j},z_{j})|t_{j}\in\mathbb{R},z_{j}\in\mathbb{C}\}$.

We now show that, if $\phi(t)=\exp\left[i\omega_{0}t-\Phi(t)\right]$ defines a CPTP pure dephasing dynamics, and if $\Phi(t)$ is even and $\phi(-t)=\phi(t)^{\ast}$, then $\phi(t)$ defined on $\mathbb{R}$ is positive definite.

Note that $\phi(t)$ describing a CPTP pure dephasing dynamics implies that $\phi(0)=1$, $\Phi(0)=0$, and $|\phi(t)|\leq\phi(0)$ for any $t>0$. This means that the coherence of the system can never exceed its initial value. These properties will be frequently used in the following proof.

Proof.

To simplify the problem, we first observe that the positive definiteness of $\phi(t)$ is equivalent to that of $\exp\left[-\Phi(t)\right]$, since

[TABLE]

Correspondingly, we can assume that $\omega_{0}=0$ without loss of generality.

Since Eq. (10) must be valid for any number of pairs, we give the proof in an inductive manner.

In the case of only one pair $(t_{1},z_{1})$, Eq. (10) is trivially satisfied. We therefore start with the case of two pairs. As stated in the main article, Eq. (10) is equivalent to the positive semidefiniteness of the Hermitian matrix:

[TABLE]

It is automatically satisfied according to the CPTP dynamics defined by $\phi(t)$.

We proceed to show the positive semidefiniteness of the Hermitian matrix

[TABLE]

for the case of three pairs. In the above matrix, and hereafter, the abbreviation $\Phi_{j,k}=\Phi(t_{j}-t_{k})$ has been adopted. Since $\mathcal{M}^{(3)}$ is three-dimensional, it is generically hard to write down an analytic expression for its three eigenvalues $\lambda_{\mu}$. Nevertheless, analyzing its characteristic polynomial gives us substantial knowledge on the eigenvalues:

1. (i)

$\lambda_{1}+\lambda_{2}+\lambda_{3}=3\geq 0$ follows from the invariance of the trace. 2. (ii)

$\lambda_{1}\lambda_{2}+\lambda_{2}\lambda_{3}+\lambda_{3}\lambda_{1}$ equals to the sum of all principal minors of $\mathcal{M}^{(3)}$ of order $2$ and is consequently non-negative, since each principal minor is non-negative, following the positive semidefiniteness of $\mathcal{M}^{(2)}$. 3. (iii)

$\lambda_{1}\lambda_{2}\lambda_{3}=\mathrm{det}\left(\mathcal{M}^{(3)}\right)$. The positivity of the product of eigenvalues is verified with the help of a simple geometric visualization shown in Fig. 1. Explicitly expanding the determinant leads to

[TABLE]

with the notation $\cos\theta=\exp\left[-\Phi_{3,2}\right]$, $a=\exp\left[-\Phi_{2,1}\right]$, and $b=\exp\left[-\Phi_{3,1}\right]$. This can be interpreted in terms of a triangle with circumcircle (dashed circle) of diameter $2r$ less than $1$. With the help of $c/\sin\theta=2r$, the positivity of Eq. (14) and, consequently, of the product of eigenvalues is then inferred.

Combining (i)-(iii), we can conclude that the three eigenvalues are non-negative each and, therefore, that $\mathcal{M}^{(3)}$ is positive semidefinite.

Before proceeding to the case of four pairs, it is worthwhile to discuss how the minimum of Eq. (10) is achieved. For the case of three pairs, the LHS of Eq. (10) is equivalent to the summation over entries in the following array:

[TABLE]

If we first determine the amplitudes $|z_{j}|$ and adjust their arguments and $t_{j}$, it is clear that the diagonal elements in the array (15) are all positive and, to reduce the resulting summation, the possible negative contributions are given by the off-diagonal elements. If we choose three pairs such that the angles between any two $z_{j}$ within them is less than $\pi/2$, as show in Fig. 2(a), the summation of all off-diagonal elements is positive. Therefore, we must choose appropriate pairs such that all their relative arguments strictly exceed $\pi/2$, as show in Fig. 2(b). To maximize the negative contributions, we assume $t_{1}=t_{2}=t_{3}$ and $\exp\left[-\Phi_{j,k}\right]=1$. We therefore draw the conclusion that $\sum_{j,k}f(t_{j}-t_{k})z_{j}z_{k}^{\ast}\geq|z_{1}+z_{2}+z_{3}|^{2}$ for the case of maximized relative arguments between three $z_{j}$.

However, in the case of four pairs, it is impossible to insert the fourth $z_{4}$ such that all relative arguments are strictly larger than $\pi/2$, as shown in Fig. 2(c). According to the above discussion, to deal with the array of four pairs,

[TABLE]

we can at most group three $z_{j}$ with all three relative arguments strictly larger than $\pi/2$ by setting their corresponding $t_{j}$ equal. Then the array (16) reduces to a simpler one:

[TABLE]

Again, in accordance with the positive semidefiniteness of $\mathcal{M}^{(2)}$, we can guarantee the validity of Eq. (10) for the case of four pairs.

For the case of five or more pairs, a similar procedure can be applied to continuously reduce the problem to an equivalent $\mathcal{M}^{(2)}$ or $\mathcal{M}^{(3)}$ case. This implies the validity of Eq. (10) for the general case. ∎

Let us remark that the above proof already indicates the general impossibility of a Hamiltonian ensemble description for arbitrary pure dephasing dynamics. Many conclusions in the above proof hold since the phase angle of $\phi(t)$ is directly proportional to time $t$. This is particularly manifest in Eq. (11). However, this is in general not the case, e.g., in the extended spin-boson model below. We consequently may obtain invalid (or quasi-) distributions in the extended spin-boson model.

III EXTENDED SPIN-BOSON MODEL

We proceed with the details of the extended spin-boson model, which consists of two qubits coupled to a common boson environment. The system and the interaction Hamiltonian of the conventional spin-boson model are thus replaced by

[TABLE]

Note that the two qubits do not interact with each other directly. Let us remark that, while we consider two qubits here, our treatment can straightforwardly be generalized to more than two qubits.

Transforming to the interaction picture with respect to $\widehat{H}_{\mathrm{S}}+\widehat{H}_{\mathrm{E}}$, the total system evolves according to the unitary evolution operator

[TABLE]

where $\mathcal{T}$ is the time-ordering operator, $\widehat{Z}_{\vec{k}}=\sum_{j=1,2}g_{j,\vec{k}}\hat{\sigma}_{z,j}$, and $\hat{b}_{\vec{k}}(t)=e^{-i\omega_{\vec{k}}t}\hat{b}_{\vec{k}}$, respectively. In the conventional spin-boson model with a single qubit, time-ordering $\mathcal{T}$ plays no significant role, since it merely introduces a global phase to the unitary evolution operator. However, this is not the case for extended models with more than one qubit, where one must carefully deal with the effect of time-ordering $\mathcal{T}$. We therefore have

[TABLE]

with

[TABLE]

By using the prescription $e^{\beta\hat{b}^{\dagger}}\hat{b}e^{-\beta\hat{b}^{\dagger}}=\hat{b}-\beta$, the operator $\widehat{\mathcal{A}}(t)$ can be recast into

[TABLE]

with

[TABLE]

Given that both $\widehat{A}$ and $\widehat{B}$ commute with $\left[\widehat{A},\widehat{B}\right]$, they satisfy $e^{\widehat{A}}e^{\widehat{B}}=e^{\left[\widehat{A},\widehat{B}\right]/2}e^{\widehat{A}+\widehat{B}}$. Then $\widehat{U}^{\mathrm{I}}(t)$ can easily be calculated:

[TABLE]

where $\alpha_{\vec{k}}(t)=-i\int_{0}^{t}e^{i\omega_{\vec{k}}\tau}d\tau=\left(1-e^{i\omega_{\vec{k}}t}\right)/\omega_{\vec{k}}$.

Assuming the direct-product initial state

[TABLE]

the reduced dynamics of qubit-1, which we now consider to be our system, can be obtained by

[TABLE]

The superscript $\mathrm{I}$ reminds that the dynamics is formulated in the interaction picture. One can easily show that the reduced dynamics of each qubit describes pure dephasing. We can thus apply the same method for constructing a Hamiltonian ensemble as for the conventional spin-boson model. We therefore focus on the time evolution of the off-diagonal element of qubit-1, which is written as

[TABLE]

where $\rho_{1,\downarrow\uparrow}(0)$, $\rho_{2,\uparrow\uparrow}(0)$, and $\rho_{2,\downarrow\downarrow}(0)$ are the initial conditions for the two qubits and the dephasing factor $\phi^{(\mathrm{X})}(t)$ is written as

[TABLE]

where

[TABLE]

$\mathcal{J}_{j,j^{\prime}}(\omega)=\sum_{\vec{k}}g_{j,\vec{k}}g_{j^{\prime},\vec{k}}^{*}\delta(\omega-\omega_{\vec{k}})$ are the spectral density functions, and

[TABLE]

represent the displacement operators, respectively.

The coupling constants $g_{j,\vec{k}}$ of the two qubits to the boson environment are in general complex numbers. In order to reveal the nonclassical effects caused by their relative phase, we assume, for simplicity, that they have the same amplitude, but with a phase difference:

[TABLE]

For a thermalized environment at temperature $T$, the two prescriptions $\exp\left(\alpha\hat{b}^{\dagger}-\alpha^{*}\hat{b}\right)\exp\left(\beta\hat{b}^{\dagger}-\beta^{*}\hat{b}\right)=\exp\left[\left(\alpha\beta^{*}-\alpha^{*}\beta\right)/2\right]\exp\left[(\alpha+\beta)\hat{b}^{\dagger}-(\alpha+\beta)^{*}\hat{b}\right]$ and $\langle\exp\left(\alpha\hat{b}^{\dagger}-\alpha^{*}\hat{b}\right)\rangle=\exp\left[-\coth(\hbar\omega/2k_{\mathrm{B}}T)|\alpha|^{2}/2\right]$, are helpful for calculating the desired result

[TABLE]

where

[TABLE]

$\mathcal{J}(\omega)=\sum_{\vec{k}}|g_{j,\vec{k}}|^{2}\delta(\omega-\omega_{\vec{k}})$ is the spectral density function, and

[TABLE]

is the same as the one in the conventional spin-boson model. In the second line of Eq. (III), we have manually inserted $\mathrm{sign}(t)$. While this does not affect the pure dephasing dynamics for $t\geq 0$, it ensures that the condition $\phi^{(\mathrm{X})}(-t)=\phi^{(\mathrm{X})*}(t)$ is satisfied and one always obtains a real distribution $\wp^{(\mathrm{X})}(\omega)$.

The presence of $\vartheta_{\varphi}(t)$ in Eq. (32) will in general result in the violation of positivity. Note that, similar to the conventional spin-boson model, individual member Hamiltonians in the Hamiltonian ensemble must be of the form $\omega\hat{\sigma}_{\mathrm{z}}/2$, which allows us to follow the same line of argument.

IV OHMIC SPECTRAL DENSITY

To demonstrate the violation of positivity explicitly, we consider the Ohmic spectral density function

[TABLE]

and the zero temperature limit where $T\rightarrow 0$. For simplicity, we also assume a degenerate system Hamiltonian.

In the case of conventional spin-boson model, the dephasing factor is

[TABLE]

and the corresponding distribution is

[TABLE]

This is obviously a legitimate probability distribution without negative values. The results are shown in Fig. 2(a) of the main article. Consequently, the Hamiltonian ensemble $\left\{\left(\omega\hat{\sigma}_{\mathrm{z}}/2,\wp_{\mathrm{o1}}(\omega)\right)\right\}$ resembles the same pure dephasing dynamics of the conventional spin-boson model characterized by $\phi_{\mathrm{o1}}(t)$. As expected, $\wp_{\mathrm{o1}}(\omega)$ is, due to the degeneracy of the system Hamiltonian, centered at $\omega=0$, and broadens with increasing $\omega_{\mathrm{c}}$.

Whereas, for the extended model, the dephasing factor reads

[TABLE]

Since the condition $\phi_{\mathrm{o1}}^{(\mathrm{X})}(-t)=\phi_{\mathrm{o1}}^{(\mathrm{X})*}(t)$ is fulfilled, the corresponding distribution $\wp_{\mathrm{o1}}^{(\mathrm{X})}(t)$ is real. However, the positivity of $\wp_{\mathrm{o1}}^{(\mathrm{X})}(t)$ is in general lost due to the presence of the nontrivial phase angle $\vartheta_{\varphi}(t)$. The results are show in Fig. 2(b) and (c) of the main article.