**1. Introduction**

In the past few decades, the discovery of fluctuation theorems [1–4] and the establishment of the framework of stochastic thermodynamics [5–7] deepened our understanding of the fluctuating nature of thermodynamic quantities (such as work, heat and entropy production) in microscopic systems [8–13]. Among various fluctuation theorems, the non-equilibrium work relation [2] sharpens our understanding of the second law of thermodynamics by presenting an elegant and precise equality associating the free energy change with the fluctuating work. Such a relation was later extended to the quantum realm based on the two-point measurement definition of the quantum fluctuating work [14,15], soon after its discovery in the classical regime. The work statistics has been widely studied in various microscopic classical and quantum systems [16–26]. Historically, the quantum– classical correspondence principle played an essential role in the development of the theory of quantum mechanics and the interpretation of the transition from quantum to classical world [27,28]. In Refs. [19,22,24], it is demonstrated that the existence of the quantum–classical correspondence principle for work distribution brings justification for the definition of quantum fluctuating work via two-point measurements.

Compared to work statistics, heat statistics relevant to thermal transport associated with a nonequilibrium stationary state has been extensively studied [29–38], but the heat statistics in a finite-time quantum thermodynamic process [39–41] and its quantum– classical correspondence have been less explored. A challenge is that the precise description of the bath dynamics requires handling a huge number of degrees of freedom of the heat bath. Different approaches have been proposed to calculate the quantum fluctuating heat and its statistics, such as the non-equilibrium Green's function approach

**Citation:** Chen, J.-F.; Qiu, T.; Quan, H.T. Quantum–Classical Correspondence Principle for Heat Distribution in Quantum Brownian Motion. *Entropy* **2021**, *23*, 1602. https://doi.org/10.3390/e23121602

Academic Editor: Ignazio Licata

Received: 30 October 2021 Accepted: 23 November 2021 Published: 29 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> School of Physics, Peking University, Beijing 100871, China; chenjinfu@pku.edu.cn (J.-F.C.); tianqiu2016@pku.edu.cn (T.Q.)

to quantum thermal transport [29,32,36,42–44] and the path-integral approach to quantum thermodynamics [45–49]. However, very few analytical results about the heat statistics have been obtained for the relaxation processes in open quantum systems. These analytical results are limited to either the relaxation dynamics described by the Lindblad master equation [39,40] or the long-time limit independent of the relaxation dynamics [50]. On the other hand, some results about the heat statistics in the classical Brownian motion model have been reported [51–59]. How the quantum and the classical heat statistics (especially associated with the relaxation dynamics in finite time) are related to each other has not been explored so far, probably due to the difficulty in studying the heat statistics in open quantum systems [60–62].

In this article, we study the heat statistics of a quantum Brownian motion model described by the Caldeira–Leggett Hamiltonian [48,63–68], where the heat bath is modeled as a collection of harmonic oscillators. Although it is well known that the dynamics of such an open quantum system can approach that of the classical Brownian motion in the classical limit *h*¯ → 0 [64], less is known about the heat statistics of this model during the finite-time relaxation process. Here, we focus on the relaxation process without external driving (the Hamiltonian of the system is time-independent); thus, the quantum fluctuating heat can be defined as the difference of the system energy between the initial and the final measurements [69].Under the Ohmic spectral density, the dynamics of the composite system is exactly solvable in the continuum limit of the bath oscillators [70]. By employing the phase-space formulation approach [71–73], we obtain analytical results of the characteristic function of heat for the Caldeira–Leggett model at any relaxation time *τ* with an arbitrary friction coefficient *κ*. Previously, such an approach was employed to study the quantum corrections to work [74–76] and entropy [77,78]. Analytical results of the heat statistics bring important insights to understand the fluctuating property of heat. By taking the classical limit *h*¯ → 0, the heat statistics of the Caldeira–Leggett model approaches that of the classical Brownian motion model. Thus, our results verify the quantum–classical correspondence principle for heat distribution, and provide justification for the definition of the quantum fluctuating heat via two-point measurements. We also verify, from the analytical results, that the heat statistics satisfies the exchange fluctuation theorem of heat [4].

The rest of this article is organized as follows. In Section 2, we introduce the Caldeira– Leggett model and define the quantum fluctuating heat. In Section 3, the analytical results of the characteristic function of heat are obtained by employing the phase-space formulation approach. We show the quantum–classical correspondence of the heat distribution and discuss the heat distribution in the long-time limit or with the extremely weak or strong coupling strength. The conclusion is given in Section 4.

#### **2. The Caldeira–Leggett Model and the Heat Statistics**

#### *2.1. The Caldeira–Leggett Model*

The quantum Brownian motion is generally described by the Caldeira–Leggett model [64,65], where the system is modeled as a single particle moving in a specific potential and the heat bath is a collection of harmonic oscillators. For simplicity, we choose the harmonic potential for the system [66,79–81], where the dynamics of such an open quantum system can be solved analytically. The system relaxes to the equilibrium state at the temperature of the heat bath. We study the heat distribution of such a quantum relaxation process and analytically obtain the characteristic function of heat and its classical correspondencebasedonthephase-spaceformulationofquantummechanics.

The total Hamiltonian of the composite system is *H*tot = *HS* + *HB* + *HSB* with each term being

$$H\_S = \frac{1}{2} \frac{\not p\_0^2}{m\_0} + \frac{1}{2} m\_0 \omega\_0^2 \hat{q}\_0^2 \tag{1}$$

$$H\_B = \sum\_{n=1}^{N} \left(\frac{1}{2}\frac{\not p\_n^2}{m\_n} + \frac{1}{2}m\_n\omega\_n^2\not q\_n^2\right) \tag{2}$$

$$H\_{SB} = -\mathfrak{q}\_0 \sum\_{n=1}^{N} \left( \mathbb{C}\_n \mathfrak{q}\_n \right) + \sum\_{n=1}^{N} \left( \frac{\mathbb{C}\_n^2}{2m\_n \omega\_n^2} \mathfrak{q}\_0^2 \right), \tag{3}$$

where *m*0, *ω*0, *q*ˆ0 and *p*ˆ0 (*mn*, *ωn*, *q*ˆ*n* and *p*ˆ*n* with *n* = 1, 2, 3, ..., *N*) are the mass, frequency, position and momentum of the system (the *n*-th bath harmonic oscillator) and *Cn* is the coupling strength between the system and the *n*-th bath harmonic oscillator. The counterterm ∑*n*[*C*<sup>2</sup> *n*/(<sup>2</sup>*mnω*<sup>2</sup> *n*)]*q*<sup>ˆ</sup> 2 0 is included in the interaction Hamiltonian *HSB* to cancel the frequency shift of the system.

The spectral density is defined as *J*(*ω*) := ∑*n*[*C*<sup>2</sup> *n*/(<sup>2</sup>*mnωn*)]*δ*(*<sup>ω</sup>* − *<sup>ω</sup>n*). We adopt an Ohmic spectral density with the Lorentz–Drude cutoff [67]

$$J(\omega) = \frac{m\_0 \kappa}{\pi} \omega \frac{\Omega\_0^2}{\Omega\_0^2 + \omega^2} \,\mathrm{}\,\tag{4}$$

where *κ* is the friction coefficient. A sufficiently large cutoff frequency Ω0 (Ω0 *ω*0) is applied to ensure a finite counter-term and the dynamics with the timescale exceeding 1/Ω0 is Markovian. Under such a spectral density, the dissipation dynamics of the Caldeira– Leggett model with a weak coupling strength *κ ω*0 reproduces that of the classical underdamped Brownian motion when taking the classical limit ¯*h* → 0 [64].

We assume the initial state to be a product state of the system and the heat bath

$$
\rho(0) = \rho\_S(0) \otimes \rho\_{\mathbb{B}'}^{\mathbb{G}}.\tag{5}
$$

which makes it possible to define the quantum fluctuating heat via two-point measurements. Here, *ρS*(0) is the initial state of the system and *ρG B* = exp(−*βHB*)/*ZB*(*β*) is the Gibbs distribution of the heat bath with the inverse temperature *β* and the partition function *ZB*(*β*) = Tr[exp(−*βHB*)].

#### *2.2. The Quantum Fluctuating Heat in the Relaxation Process*

We study the heat distribution of the relaxation process based on the two-point measurement definition of the quantum fluctuating heat. When no external driving is applied to the system, the Hamiltonian of the system is time-independent. Since no work is performed during the relaxation process, the quantum fluctuating heat can be defined as

$$Q\_{l'l} = E\_{l'}^S - E\_{l'}^S \tag{6}$$

where *E<sup>S</sup> l* (*E<sup>S</sup> l*) is the eigenenergy of the system corresponding to the outcome *l* (*l* ) at the initial (final) time *t* = 0 (*t* = *τ*). The two-point measurements over the heat bath can be hardly realized due to a huge number of degrees of freedom of the heat bath [20], while the measurements over the small quantum system are much easier in principle. The positive sign corresponds to the energy flowing from the heat bath to the system.

For the system prepared in an equilibrium state, no coherence exists in the initial state and the initial density matrix of the system commutes with the Hamiltonian of the system, [*ρ*(0), *HS*] = 0. The probability of observing the transition from *l* and *l* is

$$p\_{\pi,l'l} = \gamma\_{\pi,l'l} p\_{l'} \tag{7}$$

with the conditional transition probability *γτ*,*l<sup>l</sup>* = Tr(*P*<sup>ˆ</sup>*Sl* ⊗ *IB*)*<sup>U</sup>*tot(*τ*)(*P*<sup>ˆ</sup>*Sl* ⊗ *ρGB* )*U*†tot(*τ*) and the initial probability *pl* = Tr[*ρ*(0)*P*<sup>ˆ</sup>*Sl* ]. Here, *P*ˆ*Sl* = |*ll*| is the projection operator corresponding to the outcome *l*. The heat distribution is defined as

$$P\_{\tau}(q) := \sum\_{l',l} \delta(q - Q\_{l'l}) p\_{\tau, l'l}.\tag{8}$$

The characteristic function of heat *χτ*(*ν*) is defined as the Fourier transform of the heat distribution *χτ*(*ν*) := ∑*l*,*<sup>l</sup>* exp[*<sup>i</sup>ν*(*ESl* − *ESl* )]*p<sup>τ</sup>*,*ll*, which can be rewritten explicitly as

$$\chi\_{\tau}(\nu) = \text{Tr}\left[e^{i\nu H\_{\text{S}}} \mathcal{U}\_{\text{tot}}(\tau) \left(e^{-i\nu H\_{\text{S}}} \rho(0)\right) \mathcal{U}\_{\text{tot}}^{\dagger}(\tau)\right],\tag{9}$$

where *<sup>U</sup>*tot(*τ*) = exp(−*iH*tot*<sup>τ</sup>*/¯*h*) is the unitary time-evolution operator of the composite system.

Our goal is to analytically calculate the characteristic function *χτ*(*ν*). Previously, the quantum–classical correspondence principle for heat statistics has been analyzed with the path-integral approach to quantum thermodynamics [48], ye<sup>t</sup> the explicit result of the characteristic function (or generating function) of heat has not been obtained so far. We employ the phase-space formulation approach to solve this problem and rewrite the characteristic function Equation (9) into

$$\chi\_{\mathbf{7}}(\nu) = \text{Tr}\left[\boldsymbol{\varrho}^{i\nu H^{\rm H}\_{\rm S}(\tau)}\boldsymbol{\eta}(0)\right],\tag{10}$$

where the system Hamiltonian in the Heisenberg picture is

$$H\_S^{\rm H}(\tau) = \mathcal{U}\_{\rm tot}^{\dagger}(\tau) H\_S \mathcal{U}\_{\rm hot}(\tau),\tag{11}$$

and the density matrix-like operator *η*(0) is

$$\eta(0) = \left[ e^{-i\nu H\_S} \rho\_S(0) \right] \otimes \rho\_B^G. \tag{12}$$

We express Equation (10) with the phase-space formulation of quantum mechanics [71–76]:

$$\chi\_{\tau}(\nu) = \frac{1}{(2\pi\hbar)^{N+1}} \int d\mathbf{z} \left[ e^{i\nu H\_{\mathcal{S}}^{\rm H}(\tau)} \right]\_{\text{av}} (\mathbf{z}) \cdot P(\mathbf{z}),\tag{13}$$

where **z** represents a point **z** = [**q**, **<sup>p</sup>**]=[*q*0, ..., *qN*, *p*0, ..., *pN*] in the phase space of the composite system and the integral is performed over the whole phase space. The subscript "*w*" indicates the Weyl symbol of the corresponding operator and *<sup>P</sup>*(**z**) is the Weyl symbol of the operator *η*(0), which is explicitly defined as [71]

$$P(\mathbf{z}) := \int d\mathbf{y} \langle \mathbf{q} - \frac{\mathbf{y}}{2} \Big| \eta(0) \Big| \mathbf{q} + \frac{\mathbf{y}}{2} \rangle e^{\frac{i\mathbf{p} \cdot \mathbf{y}}{\hbar}}.\tag{14}$$

In the following, we calculate the heat statistics Equation (13) by employing the phase-space formulation approach.

#### **3. Results of the Characteristic Function of Heat**

We show a sketch of the derivation of the heat statistics *χτ*(*ν*) with the details left in Appendix A. We specifically consider the system is initially prepared at an equilibrium state *ρS*(0) = exp(−*βHS*)/*ZS*(*β*) with the inverse temperature *β* and the partition function *ZS*(*β*) = 1/[2 sinh(*βh*¯ *<sup>ω</sup>*0/2)]. The heat bath is at the inverse temperature *β*, which is different from *β*. In Equation (13), the two Weyl symbols *eiνH*H*<sup>S</sup>* (*τ*)*w*(**z**) and *<sup>P</sup>*(**z**) are obtained as

$$\left[\boldsymbol{\varepsilon}^{\dot{\boldsymbol{\mu}}\boldsymbol{\nu}\boldsymbol{H}\_{3}^{\rm th}(\tau)}\right]\_{\rm uv}(\mathbf{z}) = \frac{1}{\cos\left(\frac{\boldsymbol{v}\hbar\omega\boldsymbol{\eta}}{2}\right)} \exp\left[\frac{i}{2\hbar}\mathbf{z}^{\rm T}\mathbf{\tilde{A}}\_{\rm VZ}(\tau)\mathbf{z}\right],\tag{15}$$

and

$$P(\mathbf{z}) = \frac{2\sinh\left(\frac{\beta'\hbar\omega\_0}{2}\right)}{\cosh\left[\frac{(\beta'+i\nu)\hbar\omega\_0}{2}\right]} \cdot \left[\prod\_{n=1}^{N} 2\tanh\left(\frac{\beta\hbar\omega\_n}{2}\right)\right] \cdot \exp\left(-\frac{1}{2\hbar}\mathbf{z}^T\mathbf{A}\_{\beta\bar{z}}\mathbf{z}\right),\tag{16}$$

where the explicit expressions of the matrices **Λ ˜** *νz*(*τ*) and **<sup>Λ</sup>***βz* are given in Equations (A10) and (A36), respectively.

Substituting Equations (15) and (16) into Equation (13), the characteristic function of heat at any relaxation time *τ* with an arbitrary friction coefficient *κ* is finally obtained as

$$\chi\_{\tau}(\nu) = \left\{ \left[ (1 + i\Xi)(1 - i\Theta\Xi) - i\Xi(1 - \Theta - i\Theta\Xi) \frac{\kappa^2 \cos(2\hat{\omega}\_0 \tau) - 4\omega\_0^2}{(\kappa^2 - 4\omega\_0^2)\epsilon^{\kappa\tau}} \right]^2 \right.$$

$$+ \Xi^2 (1 - \Theta - i\Theta\Xi)^2 \left[ \left( \frac{\kappa^2 \cos(2\hat{\omega}\_0 \tau) - 4\omega\_0^2}{(\kappa^2 - 4\omega\_0^2)\epsilon^{\kappa\tau}} \right)^2 - \epsilon^{-2\kappa\tau} \right] \right\}^{\frac{1}{2}},\tag{17}$$

where the quantities Ξ and Θ are

$$\Sigma = \frac{\tan\left(\frac{\nu \hbar \omega\_0}{2}\right)}{\tanh\left[\frac{(\beta' + i\nu)\hbar \omega\_0}{2}\right] - i \tan\left(\frac{\nu \hbar \omega\_0}{2}\right)},\tag{18}$$

$$\Theta = \frac{\tanh\left[\frac{(\beta'+i\nu)\hbar\omega\_0}{2}\right] - i\tan\left(\frac{\nu\hbar\omega\_0}{2}\right)}{\tanh\left(\frac{\beta\hbar\omega\_0}{2}\right)}.\tag{19}$$

Induced by the friction, the frequency of the system harmonic oscillator is shifted to *ω* ˆ 0 = *ω*<sup>20</sup> − *κ*2/4.

From the analytical results of the heat statistics Equation (17), the average heat *Q*-(*τ*) = −*i∂ν*[ln *χτ*(*ν*)]|*<sup>ν</sup>*=<sup>0</sup> is immediately obtained as

$$\langle Q \rangle(\tau) = \frac{\omega\_0 \hbar}{2} \left[ \coth \left( \frac{\beta \omega\_0 \hbar}{2} \right) - \coth \left( \frac{\beta' \omega\_0 \hbar}{2} \right) \right] \left[ 1 - \frac{\kappa^2 \cos(2\omega\_0 \tau) - 4\omega\_0^2}{(\kappa^2 - 4\omega\_0^2)\epsilon^{\kappa \tau}} \right], \tag{20}$$

and the variance Var(*Q*)(*τ*) = <sup>−</sup>*∂*2*ν*[ln *χτ*(*ν*)]-*<sup>ν</sup>*=0 is

$$\text{Var}(Q)(\tau) = \text{I} + \text{II} \cdot e^{-\kappa \tau} + \text{III} \cdot e^{-2\kappa \tau},\tag{21}$$

with

$$\mathbf{I} = \frac{\omega\_0^2 \hbar^2 \left[ \cosh^2 \left( \frac{\beta \omega\_0 \hbar}{2} \right) + \cosh^2 \left( \frac{\beta' \omega\_0 \hbar}{2} \right) \right]}{4},\tag{22}$$

$$\Pi = \frac{\kappa^2 \cos(2\omega\_0 \tau) - 4\omega\_0^2}{2\omega\_0^2} \cdot \frac{\omega\_0^2 \hbar^2 \left[ \coth^2\left(\frac{\beta \omega\_0 \hbar}{2}\right) + \csc\hbar^2 \left(\frac{\beta' \omega\_0 \hbar}{2}\right) - \coth\left(\frac{\beta \omega\_0 \hbar}{2}\right) \coth\left(\frac{\beta' \omega\_0 \hbar}{2}\right) \right]}{4} \tag{23}$$

$$\mathrm{III} = \frac{\kappa^4 \cos(4\hat{\omega}\_0 \tau) + 8\omega\_0^2 \kappa^2 [1 - 2\cos(2\hat{\omega}\_0 \tau)] + 16\omega\_0^4}{16\omega\_0^4} \cdot \frac{\omega\_0^2 \hbar^2 \left[\coth\left(\frac{\beta \omega\_0 \hbar}{2}\right) - \coth\left(\frac{\beta' \omega\_0 \hbar}{2}\right)\right]^2}{4}. \tag{24}$$

Similarly, one can calculate the higher cumulants from the analytical results of the heat statistics. In the following, we examine the properties of the heat statistics of the quantum Brownian motion.

#### *3.1. Quantum–Classical Correspondence Principle for Heat Statics and the Exchange Fluctuation Theorem of Heat*

We further take the classical limit *h*¯ → 0 or, more rigorously, *βh*¯ *ω*0 → 0. The two quantities approach Ξ → *ν*/*β* and Θ → *β*/*β* and the characteristic function of heat (Equation (17)) becomes

$$\chi\_{\tau}^{\rm cl}(\nu) = \left\{ \left[ (1 + i\frac{\nu}{\beta'})(1 - i\frac{\nu}{\beta}) - \frac{i\nu(\beta - \beta' - i\nu)}{\beta\beta'} \frac{(\kappa^2 \cos(2\omega\_0 \tau) - 4\omega\_0^2)}{(\kappa^2 - 4\omega\_0^2)\epsilon^{\kappa\tau}} \right]^2 \right.$$

$$+ \nu^2 \left( \frac{\beta - \beta' - i\nu}{\beta\beta'} \right)^2 \left[ \left( \frac{\kappa^2 \cos(2\omega\_0 \tau) - 4\omega\_0^2}{(\kappa^2 - 4\omega\_0^2)\epsilon^{\kappa\tau}} \right)^2 - \epsilon^{-2\kappa\tau} \right] \right\}^{-\frac{1}{2}},\tag{25}$$

which is consistent with the results obtained from the classical Brownian motion described by the Kramers equation (see Ref. [58] or Appendix C). The average heat is

$$
\left< \mathcal{Q}^{\rm cl} \right> (\tau) = \frac{\beta' - \beta}{\beta \beta'} \left[ 1 - \frac{\kappa^2 \cos(2\omega\_0 \tau) - 4\omega\_0^2}{(\kappa^2 - 4\omega\_0^2)e^{\kappa \tau}} \right], \tag{26}
$$

and the variance Var*Q*cl(*τ*) = <sup>−</sup>*∂*2*ν*[ln *χ*cl *τ* (*ν*)]--*<sup>ν</sup>*=0 is

$$\text{Var}\left(\boldsymbol{Q}^{\text{cl}}\right)(\boldsymbol{\tau}) = \mathbf{I}^{\text{cl}} + \boldsymbol{\Pi}^{\text{cl}} \cdot \boldsymbol{\varepsilon}^{-\kappa \tau} + \mathbf{III}^{\text{cl}} \cdot \boldsymbol{\varepsilon}^{-2\kappa \tau},\tag{27}$$

with

$$
\Lambda^{\rm cl} = \frac{\beta^2 + \beta'^2}{\beta^2 \beta'^2} \tag{28}
$$

$$\mathrm{II}^{\mathrm{cl}} = \frac{\kappa^2 \cos(2\omega\_0 \tau) - 4\omega\_0^2}{2\omega\_0^2} \cdot \frac{\beta^2 - \beta \beta' + \beta'^2}{\beta^2 \beta'^2} \tag{29}$$

$$\mathrm{III}^{\mathrm{cl}} = \frac{\kappa^4 \cos(4\hat{\omega}\_0 \tau) + 8\omega\_0^2 \kappa^2 [1 - 2\cos(2\hat{\omega}\_0 \tau)] + 16\omega\_0^4}{16\omega\_0^4} \cdot \frac{(\beta - \beta')^2}{\beta^2 \beta'^2} . \tag{30}$$

 From Equation (17) (or the classical counterpart Equation (25)), one can see the characteristic function of heat exhibits the following symmetry:

$$\chi\_{\mp}(\nu) = \chi\_{\mp}[-i(\beta - \beta') - \nu],\tag{31}$$

which shows that the heat distribution satisfies the exchange fluctuation theorem of heat in the differential form *<sup>P</sup>τ*(*Q*)/*Pτ*(−*Q*) = exp[−(*β* − *β*)*Q*] [4]. By setting *ν* = 0, we obtain the relation *χτ*[−*<sup>i</sup>*(*β* − *β*)] = *χτ*(0) = 1, which is exactly the exchange fluctuation theorem of heat in the integral form exp[(*β* − *β*)*Q*]- = 1.

### *3.2. Long-Time Limit*

In the long-time limit *τ* → <sup>∞</sup>, the characteristic functions of heat (Equations (17) and (25)) become

$$\chi\_{\infty}(\nu) = \frac{\left(1 - e^{-\beta'\omega\_0\hbar}\right)\left(1 - e^{-\beta\omega\_0\hbar}\right)}{\left(1 - e^{-\left(\beta' + i\nu\right)\omega\_0\hbar}\right)\left(1 - e^{-\left(\beta - i\nu\right)\omega\_0\hbar}\right)},\tag{32}$$

and

$$\chi^{\rm cl}\_{\infty}(\nu) = \frac{\beta^{\prime}\beta}{(\beta^{\prime} + i\nu)(\beta - i\nu)}.\tag{33}$$

Such results, independent of the relaxation dynamics, are in the form

$$\chi\_{\rm th}(\nu) = \frac{Z\_{\mathcal{S}}(\boldsymbol{\beta}' + i\nu)Z\_{\mathcal{S}}(\boldsymbol{\beta} - i\nu)}{Z\_{\mathcal{S}}(\boldsymbol{\beta}')Z\_{\mathcal{S}}(\boldsymbol{\beta})},\tag{34}$$

reflecting complete thermalization of the system [53]. For example, the relaxation of a harmonic oscillator governed by the quantum–optical master equation gives the identical characteristic function of heat in the long-time limit [39]. In Appendix D, we demonstrate that the characteristic function of heat for any relaxation process with complete thermalization is always in the form of Equation (34). With the simple expressions (32) and (33) of the characteristic functions, the heat distributions are obtained from the inverse Fourier transform as

$$P\_{\infty}(q) = \begin{cases} \frac{\left(1 - e^{-\beta'\omega\_0\hbar}\right)\left(1 - e^{-\beta\omega\_0\hbar}\right)}{1 - e^{-\left(\beta'+\beta\right)\omega\_0\hbar}} \sum\_{j=0}^{\infty} \delta\left(q - j\omega\_0\hbar\right) e^{-\beta q} & q \ge 0\\ \frac{\left(1 - e^{-\beta'\omega\_0\hbar}\right)\left(1 - e^{-\beta\omega\_0\hbar}\right)}{1 - e^{-\left(\beta'+\beta\right)\omega\_0\hbar}} \sum\_{j=1}^{\infty} \delta\left(q + j\omega\_0\hbar\right) e^{\beta'q} & q < 0 \end{cases} \tag{35}$$

and

$$P\_{\infty}^{\rm cl}(q) = \begin{cases} \frac{\beta^{\prime}\beta}{\beta^{\prime} + \beta^{\prime}} e^{-\beta q} & q \ge 0\\ \frac{\beta^{\prime}\beta}{\beta^{\prime} + \beta^{\prime}} e^{\beta^{\prime}q} & q < 0 \end{cases} \tag{36}$$

which are exactly the same as the long-time results obtained in Ref. [39].

#### *3.3. Weak/Strong Coupling Limit in Finite Time*

In the weak coupling limit *κ ω*0, the characteristic function of heat Equation (17) becomes

$$\chi^{\rm w}\_{\tau}(\nu) = \frac{1}{(1 + i\Xi)(1 - i\Xi\Theta)(1 - e^{-\kappa\tau}) + e^{-\kappa\tau}}.\tag{37}$$

There is only one relaxation timescale associated to *κ*. Such situation corresponds to the highly underdamped regime of the classical Brownian motion and a systematic method has been proposed to study the heat distribution [56], as well as the work distribution, under an external driving [82,83].

In the strong coupling limit *κ ω*0, the characteristic function of heat Equation (17) becomes

$$\begin{split} \chi^{s}\_{\tau}(\nu) &= \frac{1}{\sqrt{(1+i\Xi)(1-i\Xi\Theta)(1-e^{-2\chi\tau}) + e^{-2\chi\tau}}} \\ &\times \frac{1}{\sqrt{(1+i\Xi)(1-i\Xi\Theta)\left(1-e^{-\frac{2\omega\_{0}^{2}}{\kappa}\tau}\right) + e^{-\frac{2\omega\_{0}^{2}}{\kappa}\tau}}}. \tag{38} \end{split} \tag{38}$$

The relaxation timescales of the momentum (the first factor) and the coordinate (the second factor) are separated. The long-time limits of both Equations (37) and (38) are equal to Equation (32). In classical thermodynamics, the usual overdamped approximation neglects the motion of the momentum; hence, the heat statistics derived under such an approximation is incomplete [52]. Actually, the momentum degree of freedom also contributes to the heat statistics.

#### *3.4. Numerical Results*

In Figure 1, we show the cumulative heat distribution function Pr(*Q* < *q*) := *q* −∞ *<sup>P</sup>τ*(*q* )*dq* with different friction coefficients *κ* = 0.01, 1 and 100, at the rescaled relaxation time *τ*˜ = *κτ* = 1 and 10. We set the mass *m*0 = 1 and the frequency *ω*0 = 1 for the system harmonic oscillator, the inverse temperatures *β* = 1 and *β* = 2 for the initial equilibrium states of the heat bath and the system, respectively. The Planck constant is set to be *h*¯ = 1, 0.5, 0.1. With the decrease in *h*¯, the quantum result Equation (17) approaches the classical result Equation (25). Thus, the quantum–classical correspondence of the heat distribution is demonstrated for generic values of the friction coefficient *κ*.

**Figure 1.** The cumulative heat distribution function Pr(*Q* < *q*). The choices of the parameters are given in the main text. We compare the results of the Caldeira–Leggett model (blue solid, orange dotted and green dot-dashed curves) in Equation (17) and those of the classical Brownian motion (black dashed curve) in Equation (25). The rescaled relaxation time is *τ*˜ = *κτ* = 1 in the upper subfigures and *τ*˜ = 10 in the lower subfigures. The left, middle and right subfigures illustrate the results for the weak (*κ* = 0.01), intermediate (*κ* = 1) and strong coupling strength (*κ* = 100).

For *κ* = 0.01 and 1, complete thermalization is achieved at *τ*˜ = 10. The leftlower and middle-lower subfigures show the identical distribution characterized by Equations (35) and (36). For *κ* = 100, the momentum degree of freedom is thermalized exp(−2*τ*˜) ≈ 0 in Equation (38), while the coordinate degree of freedom remains frozen exp[−<sup>2</sup>(*ω*<sup>2</sup> 0/*κ*<sup>2</sup>)*τ*˜] ≈ 1 in Equation (38). Thus, the distribution in the right-lower subfigure is different from the middle-lower subfigure.

In Figure 2, we illustrate the results of the mean value *Q*-(*τ*) and the variance Var(*Q*)(*τ*) with different friction coefficients *κ* = 0.01, 1 and 100. The parameters are the same as those in Figure 1. The quantum results approach the classical results with the decrease in *h*¯. For *κ* = 0.01 and 1 (left and middle subfigures), complete thermalization is reached when *τ*˜ > 5. The mean value and the variance approach lim*τ*→ ∞ *Q*cl (*τ*) = 1/*β* − 1/*β* and lim*τ*→ ∞ Var *Q*cl (*τ*) = 1/*β*<sup>2</sup> + 1/*β*2(gray horizontal lines). For *κ* = 100 (right subfigures), only the momentum degree of freedom is thermalized at this timescale. Thus, the mean value and the variance take half value of their long-time limits. When the

coordinate degree of freedom is also thermalized in the long-time limit (*τ*˜ *κ*2/*ω*20 = 104), the mean value and the variance are expected to approach the same values as those in the middle subfigures.

**Figure 2.** The evolution of the mean value *Q*-(*τ*) (**upper subfigures**) and the variance Var(*Q*)(*τ*) (**lower subfigures**) of the heat statistics as functions of the rescaled time *τ*˜ = *κτ*.

### **4. Conclusions**

Previously, the heat statistics of the relaxation processes has been studied analytically in open quantum systems described by the Lindblad master equation [39,40,50]. However, due to the rotating wave approximation and other approximations, such quantum systems do not possess a well-defined classical counterpart. Hence, the quantum–classical correspondence principle for heat distribution has not been well established.

In this paper, we study the heat statistics of the quantum Brownian motion model described by the Caldeira–Leggett Hamiltonian, in which the bath dynamics is explicitly considered. By employing the phase-space formulation approach, we obtain the analytical expressions of the characteristic function of heat at any relaxation time *τ* with an arbitrary friction coefficient *κ*. The analytical results of heat statistics bring important insights to the studies of quantum thermodynamics. For example, in the classical limit, our results approach the heat statistics of the classical Brownian motion. Thus, the quantum–classical correspondence principle for heat statistics is verified in this model. Our analytical results provide justification for the definition of quantum fluctuating heat via two-point measurements.

We also discuss the characteristic function of heat in the long-time limit or with the extremely weak/strong coupling strength. In the long-time limit, the form of the characteristic function of heat reflects complete thermalization of the system. In addition, from the analytical expressions of the heat statistics, we can immediately verify the exchange fluctuation theorem of heat. The phase-space formulation can be further utilized to study the joint statistics of work and heat in a driven open quantum system, which would be beneficial in exploring the fluctuations of power and efficiency in finite-time quantum heat engines.

**Author Contributions:** Conceptualization, J.-F.C. and H.-T.Q.; Formal analysis, J.-F.C.; Funding acquisition, H.-T.Q.; Investigation, J.-F.C.; Methodology, J.-F.C.; Supervision, H.-T.Q.; Validation, T.Q. and H.-T.Q.; Visualization, J.-F.C.; Writing—original draft, J.-F.C.; Writing—review & editing, T.Q. and H.-T.Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** H. T. Quan acknowledges support from the National Natural Science Foundation of China under Grants No. 11775001, No. 11534002 and No. 11825001.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This paper is dedicated to Wojciech Zurek on the occasion of his 70th birthday for his mentorship as well as his kind and generous supports to one of the authors (H. T. Quan), and for his many insightful contributions to our understanding about the quantum-to-classical transition.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Derivation to the Characteristic Function of Heat (17)**

We show the detailed derivation to the characteristic function of heat *χτ*(*ν*). First, we calculate the two Wigner functions *eiνH*H*<sup>S</sup>* (*t*)*w*(**z**) and *<sup>P</sup>*(**z**). Then, the final result Equation (17) is obtained from Equation (13).

$$\operatorname{Appendix } \mathcal{A}. \mathbf{1.} \left[ \mathbf{e}^{i\nu H\_S^H(t)} \right]\_w (\mathbf{z})$$

With the quadratic Hamiltonian *<sup>H</sup>*H*S* (*t*), the Wigner function *eiνH*H*<sup>S</sup>* (*t*)*w*(**z**) is [78,81]

$$\left[\epsilon^{j\nu H\_3^{\rm H}(t)}\right]\_{w}(\mathbf{z}) = \frac{1}{\cos\left(\frac{\omega\_0 \hbar \nu}{2}\right)} \exp\left[i\frac{m\_0 \omega\_0}{\hbar} \tan\left(\frac{\omega\_0 \hbar \nu}{2}\right) q\_0^2(t) + i \frac{1}{m\_0 \hbar \omega\_0} \tan\left(\frac{\omega\_0 \hbar \nu}{2}\right) p\_0^2(t)\right] \tag{A1}$$

$$\mathbf{x} = \frac{1}{\cos\left(\frac{\omega\_0 \hbar v}{2}\right)} \exp\left[\frac{i}{2\hbar} \mathbf{z}^T(t) \mathbf{A}\_{vz} \mathbf{z}(t)\right],\tag{A2}$$

where **z**(*t*) gives the trajectory in the phase space determined by the initial point **z**(0) = **z** and **Λ***νz* is a rank-2 diagonal matrix

$$\mathbf{A}\_{\upsilon z} = \begin{pmatrix} 2m\_0 \omega\_0 \tan\left(\frac{\omega\_0 \hbar \upsilon}{2}\right) \\\\ \mathbf{O} \\\\ \frac{2}{m\_0 \omega\_0} \tan\left(\frac{\omega\_0 \hbar \upsilon}{2}\right) \\\\ \mathbf{O} \end{pmatrix}' \tag{A3}$$

with an *N* × *N* zero matrix **O**. The unlisted elements are zeros. The trajectory **z**(*t*) satisfies the classical equation of motion (also the equation of motion in the Heisenberg picture).

$$
\dot{q}\_0 = \frac{p\_0}{m\_0} \,\tag{A4}
$$

$$\dot{q}\_n = \frac{p\_n}{m\_n}'\tag{A5}$$

$$
\dot{p}\_0 = -m\_0 \omega\_0^2 q\_0 + \sum\_n \mathbb{C}\_n q\_n,\tag{A6}
$$

$$
\dot{p}\_n = -m\_n \omega\_n^2 q\_n + \mathcal{C}\_n q\_{0\prime} \tag{A7}
$$

with *ω*˜ 20 = *ω*20 + ∑*Nn*=<sup>1</sup> *<sup>C</sup>*2*n*/(*<sup>m</sup>*0*mnω*2*n*). The above differential equations can be rewritten into a compact form **z**˙(*t*) = **Lz**(*t*). The trajectory **<sup>z</sup>**(*t*)=[**q**(*t*), **p**(*t*)] in the phase space characterizes the evolution of the composite system with the positions **<sup>q</sup>**(*t*)=[*q*0(*t*), ..., *qN*(*t*)] and the momenta **<sup>p</sup>**(*t*)=[*p*0(*t*), ..., *pN*(*t*)] and is related to the initial point by the dynamical map **z**(*t*) = exp(**<sup>L</sup>***t*)**z**(0). The (2*N* + 2) × (2*N* + 2) matrix **L** is explicitly

**L** = ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1*m*0 1*m*1 1*m*2 ... 1*mN* −*m*0*ω*˜ 20 *C*1 *C*2 ... *CN C*1 −*m*1*<sup>ω</sup>*21 *C*2 −*m*2*<sup>ω</sup>*22 ... ... *CN* <sup>−</sup>*mNω*2*N* ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠, (A8)

and the matrix exponential is formally written as

$$\exp(\mathbf{L}t) = \begin{pmatrix} a\_0 & a\_1 & \dots & a\_N & \frac{\beta\_0}{m\_0} & \frac{\beta\_1}{m\_1} & \dots & \frac{\beta\_N}{m\_N} \\ \gamma\_1 & \Lambda\_{11} & \dots & \Lambda\_{1N} & \frac{\xi\_1}{m\_0} & \frac{\Lambda\_{11}}{m\_1} & \dots & \frac{\Lambda\_{1N}}{m\_N} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \gamma\_N & \Lambda\_{N1} & \dots & \Lambda\_{NN} & \frac{\xi\_N}{m\_0} & \frac{\beta\_{N1}}{m\_1} & \dots & \frac{\Lambda\_{NN}}{m\_N} \\ m\_0\mu\_0 & m\_0\mu\_1 & \dots & m\_0\mu\_N & \frac{m\_0}{m\_0}\beta\_0 & \frac{m\_0}{m\_1}\beta\_1 & \dots & \frac{m\_0}{m\_N}\beta\_N \\ m\_1\dot{\gamma}\_1 & m\_1\Lambda\_{11} & \dots & m\_1\Lambda\_{1N} & \frac{m\_1}{m\_0}\xi\_1 & \frac{m\_1}{m\_1}\dot{\Lambda}\_{11} & \dots & \frac{m\_1}{m\_N}\Lambda\_{1N} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ m\_N\dot{\gamma}\_N & m\_N\Lambda\_{N1} & \dots & m\_N\Lambda\_{NN} & \frac{m\_N}{m\_0}\xi\_N & \frac{m\_N}{m\_1}\Lambda\_{N1} & \dots & \frac{m\_N}{m\_N}\dot{\Lambda}\_{NN} \end{pmatrix} . \tag{A9}$$

We rewrite the quadratic form into **<sup>z</sup>**<sup>T</sup>(*t*)**<sup>Λ</sup>***νz***<sup>z</sup>**(*t*) = **z**<sup>T</sup>(0)**Λ˜** *νz*(*t*)**z**(0) with

$$\tilde{\mathbf{A}}\_{\nu z}(t) = \exp\left(\mathbf{L}^{\mathrm{T}}t\right) \mathbf{A}\_{\nu z} \exp(\mathbf{L}t). \tag{A10}$$

We next carry out every element in Equation (A9) through the Laplace transforms of Equations (A4)–(A7):

$$s\overline{q}\_0(s) - q\_0(0) = \frac{\overline{p}\_0(s)}{m\_0},\tag{A11}$$

$$\text{s}\ddot{q}\_{\text{n}}(\text{s}) - q\_{\text{n}}(\text{0}) = \frac{\tilde{p}\_{\text{n}}(\text{s})}{m\_{\text{n}}},\tag{A12}$$

$$s\overline{p}\_0(s) - p\_0(0) = -m\_0 \tilde{\omega}\_0^2 \overline{q}\_0(s) + \sum\_n \mathbb{C}\_n \overline{q}\_n(s),\tag{A13}$$

$$s\overline{p}\_n(s) - p\_n(0) = -m\_n \omega\_n^2 \overline{q}\_n(s) + \mathbb{C}\_n \overline{q}\_0(s). \tag{A14}$$

Representing *q*˜*n*(*s*) and *p*˜*n*(*s*) with *q*˜0(*s*) and the initial conditions, we obtain

$$\left\{s^{2}+\omega\_{0}^{2}-\sum\_{n}\left[\frac{\mathcal{C}\_{n}^{2}}{m\_{0}m\_{n}(s^{2}+\omega\_{n}^{2})}\right]\right\}\dot{q}\_{0}(s)=\dot{q}\_{0}(0)+sq\_{0}(0)+\sum\_{n}\frac{\mathcal{C}\_{n}}{m\_{0}}\left[\frac{\dot{q}\_{n}(0)+sq\_{n}(0)}{s^{2}+\omega\_{n}^{2}}\right].\tag{A15}$$

Under the Ohmic spectral density Equation (4), the above equation is simplified to

$$(s^2 + \kappa s + \omega\_0^2)\ddot{q}\_0(s) = \dot{q}\_0(0) + s\dot{q}\_0(0) + \sum\_n \frac{\mathbb{C}\_n}{m\_0} \left[ \frac{\dot{q}\_n(0) + sq\_n(0)}{s^2 + \omega\_n^2} \right],\tag{A16}$$

where the summation on the left-hand side of Equation (A15) can be approximately expressed as

$$\sum\_{n} \left[ \frac{\mathcal{C}\_{n}^{2}}{m\_{0}m\_{n} \left(s^{2} + \omega\_{n}^{2}\right)} \right] \approx -\kappa \mathbf{s} + \sum\_{n} \frac{\mathcal{C}\_{n}^{2}}{m\_{0}m\_{n}\omega\_{n}^{2}},\tag{A17}$$

with a large cutoff frequency Ω0. The inverse Laplace transform gives the differential equation of *q*0(*t*) as

$$\ddot{q}\_0(t) + \kappa \dot{q}\_0(t) + \omega\_0^2 q\_0(t) = \underbrace{-\kappa q\_0(0)\delta(t)}\_{\text{initial velocity change}} + \underbrace{\sum\_n \frac{\mathbb{C}\_n}{m\_0} \left[ \dot{q}\_n(0) \frac{\sin(\omega\_n t)}{\omega\_n} + q\_n(0) \cos(\omega\_n t) \right]}\_{\text{stochastic force}}.\tag{A18}$$

On the right-hand side, the second term presents the stochastic force induced by the heat bath; the first term indicates an abrupt velocity change −*<sup>κ</sup>q*0(0) of the system particle at the initial time *t* = 0 [63,84,85]. The sudden change in velocity occurs for the system harmonic oscillator when the coupling between the system and the heat bath is switched on. Such an initial slippage is caused by the assumption of the initial product state. To avoid such an initial discontinuous problem, we drop the first term by considering the particle motion as starting at *t* = 0+ [70]. Under such a modification, the Caldeira–Leggett model can reproduce the complete Langevin equation with an arbitrary friction coefficient *κ* for both the underdamped and the overdamped regimes and the heat distribution of the Caldeira–Leggett model approaches that of the classical Brownian motion described by the Kramers equation [86]. In Appendix B, for the classical counterpart of the Caldeira–Leggett model, we show that the initial slippage can be naturally eliminated by choosing another initial state.

After dropping the first term, Equation (A16) becomes

$$\dot{q}\_{1}(\mathbf{s}^{2}+\mathbf{x}\mathbf{s}+\omega\_{0}^{2})\ddot{q}\_{0}(\mathbf{s})=\dot{q}\_{0}(0)+(\mathbf{x}+\mathbf{s})q\_{0}(0)+\sum\_{n}\frac{\mathbf{C}\_{n}}{m\_{0}}\left[\frac{\dot{q}\_{n}(0)+\mathbf{s}q\_{n}(0)}{\mathbf{s}^{2}+\omega\_{n}^{2}}\right].\tag{A19}$$

The solutions to *q*˜0(*s*) and *q*˜*n*(*s*) follow immediately as

$$\vec{q}\_0(s) = \frac{\dot{q}\_0(0) + (\kappa + s)q\_0(0) + \sum\_n \frac{\mathbb{C}\_n}{m\_0} \frac{q\_n(0) + sq\_n(0)}{s^2 + \omega\_n^2}}{s^2 + \kappa s + \omega\_0^2},\tag{A20}$$

$$\vec{q}\_{n}(\mathbf{s}) = \frac{\dot{q}\_{n}(0) + sq\_{n}(0)}{\mathbf{s}^{2} + \omega\_{n}^{2}} + \frac{\mathbb{C}\_{n}}{m\_{n}} \cdot \frac{\dot{q}\_{0}(0) + (\mathbf{x} + \mathbf{s})q\_{0}(0) + \sum\_{l} \frac{\mathbb{C}\_{l}}{m\_{l}} \frac{\dot{q}\_{l}(0) + sq\_{l}(0)}{\mathbf{s}^{2} + \omega\_{l}^{2}}}{\left(\mathbf{s}^{2} + \omega\_{n}^{2}\right)\left(\mathbf{s}^{2} + \kappa\mathbf{s} + \omega\_{0}^{2}\right)}.\tag{A21}$$

With the inverse Laplace transform, the elements in the matrix exp(**<sup>L</sup>***t*) Equation (A9) are determined by

$$
\begin{pmatrix} q\_0(t) \\ q\_1(t) \\ \dots \\ q\_N(t) \end{pmatrix} = \begin{pmatrix} a\_0 & a\_1 & a\_2 & \dots & a\_N & \beta\_0 & \beta\_1 & \beta\_2 & \dots & \beta\_N \\ \gamma\_1 & \Lambda\_{11} & \Lambda\_{12} & \dots & \Lambda\_{1N} & \xi\_1 & \Delta\_{11} & \Delta\_{12} & \dots & \Delta\_{1N} \\ \dots & \dots & & \dots & & \dots & & \dots \\ \gamma\_N & \Lambda\_{N1} & \Lambda\_{N2} & \dots & \Lambda\_{NN} & \xi\_N & \Delta\_{N1} & \Delta\_{N2} & \dots & \Delta\_{NN} \end{pmatrix} \begin{pmatrix} q\_0(0) \\ q\_1(0) \\ \dots \\ q\_N(0) \\ \dot{q}\_0(0) \\ \vdots \\ \dot{q}\_1(0) \\ \vdots \\ \dot{q}\_N(0) \end{pmatrix}, \tag{A22}
$$

where the elements in the matrix of the right-hand side are explicitly solved as [70]

$$\boldsymbol{\alpha}\_{0} = \boldsymbol{e}^{-\frac{\boldsymbol{\kappa}!}{2}} \left[ \cos(\hat{\boldsymbol{\omega}}\_{0} \boldsymbol{t}) + \frac{\boldsymbol{\kappa}}{2\hat{\boldsymbol{\omega}}\_{0}} \sin(\hat{\boldsymbol{\omega}}\_{0} \boldsymbol{t}) \right],\tag{A23}$$

$$\beta\_0 = \frac{e^-\pi}{\hat{\omega}\_0} \sin(\hat{\omega}\_0 t),\tag{A24}$$

$$\alpha\_{\rm n} = \frac{\mathcal{C}\_{\rm n}}{m\_0} f\_{\rm n}(t),\tag{A25}$$

$$
\beta\_n = \frac{\mathbb{C}\_n}{m\_0} \mathbf{g}\_n(t),
\tag{A26}
$$

$$\gamma\_n = \frac{\mathbb{C}\_n}{m\_n} [f\_n(t) + \kappa \mathbb{g}\_n(t)],\tag{A27}$$

$$\mathbf{g}\_{\text{fl}}^{\text{g}} = \frac{\mathbb{C}\_{\text{n}}}{m\_{\text{n}}} \mathbf{g}\_{\text{n}}(t),\tag{A28}$$

$$
\Lambda\_{nl} = \delta\_{nl}\cos(\omega\_n t) + \frac{\mathcal{C}\_n \mathcal{C}\_l}{m\_n m\_0} F\_{nl}(t), \tag{A29}
$$

$$
\Delta\_{nl} = \frac{\delta\_{nl}}{\omega\_{ll}} \sin(\omega\_{ll} t) + \frac{\mathcal{C}\_{\text{n}} \mathcal{C}\_{l}}{m\_{\text{n}} m\_{\text{0}}} \mathcal{G}\_{\text{nl}}(t). \tag{A30}
$$

The functions *fn*(*t*), *gn*(*t*), *Fnl*(*t*) and *Gnl*(*t*) are, explicitly,

$$f\_n(t) = \mathcal{L}^{\ell-1} \left[ \frac{\mathsf{s}}{(\mathsf{s}^2 + \mathsf{k}\mathsf{s} + \omega\_0^2)(\mathsf{s}^2 + \omega\_n^2)} \right]\_{\mathsf{n}}^{\prime} \tag{A31}$$

$$g\_n(t) = \mathcal{E}^{\prime -1} \left[ \frac{1}{(s^2 + \kappa s + \omega\_0^2)(s^2 + \omega\_n^2)} \right],\tag{A32}$$

$$F\_{nl}(t) = \mathcal{L}^{\rho - 1} \left[ \frac{\mathbf{s}}{(s^2 + \kappa \mathbf{s} + \omega\_0^2)(s^2 + \omega\_n^2)(s^2 + \omega\_l^2)} \right]\_{\mathbf{n}}^{\mathbf{n}} \tag{A33}$$

$$G\_{nl}(t) = \mathcal{L}^{\rho - 1} \left[ \frac{1}{(s^2 + \kappa s + \omega\_0^2)(s^2 + \omega\_n^2)(s^2 + \omega\_l^2)} \right],\tag{A34}$$

where L <sup>−</sup><sup>1</sup>(·) denotes the inverse Laplace transform with L (·) = ∞0 (·)*e*<sup>−</sup>*stdt*.

### *Appendix A.2. <sup>P</sup>*(**z**)

*<sup>P</sup>*(**z**) is the Wigner function of the state *η*(0) for the composite system [78,81]:

$$P(\mathbf{z}) = \frac{2\sinh\left(\frac{\beta^{\prime}\hbar\omega\_{0}}{2}\right)}{\cosh\left[\frac{(\beta^{\prime}+i\nu)\hbar\omega\_{0}}{2}\right]} \cdot \left[\prod\_{n=1}^{N} 2\tanh\left(\frac{\beta\hbar\omega\_{n}}{2}\right)\right] \cdot \exp\left[-\frac{1}{2\hbar}\mathbf{z}^{\mathrm{T}}\mathbf{A}\_{\beta\mathrm{z}}\mathbf{z}\right],\tag{A35}$$

where **<sup>Λ</sup>***βz* is a (2*N* + 2) × (2*N* + 2) diagonal matrix

$$\mathbf{A}\_{\beta z} = \text{diag}\left(\lambda\_{\beta' \neq 0'} \theta\_1, \dots, \theta\_{N'} \lambda\_{\beta' p\_{0'}} \mu\_{1'}, \dots, \mu\_N\right), \tag{A36}$$

with the elements

$$\theta\_n = 2m\_n \omega\_n \tanh\left(\frac{\beta \hbar \omega\_n}{2}\right),\tag{A37}$$

$$
\mu\_{\rm ll} = \frac{2}{m\_{\rm n}\omega\_{\rm n}} \tanh\left(\frac{\beta\hbar\omega\_{\rm n}}{2}\right),
\tag{A38}
$$

$$
\lambda\_{\beta'q\_0} = 2m\_0 \omega\_0 \tanh\left[\frac{(\beta'+i\nu)\hbar\omega\_0}{2}\right],\tag{A39}
$$

$$
\lambda\_{\beta'p\_0} = \frac{2}{m\_0 \omega\_0} \tanh\left[\frac{(\beta' + i\nu)\hbar\omega\_0}{2}\right].\tag{A40}
$$

#### *Appendix A.3. Calculation of the Integral*

With the explicit expressions of *eiνH*H*<sup>S</sup>* (*τ*)*w*(**z**) and *<sup>P</sup>*(**z**), we perform the integral in Equation (13) and obtain the result of the characteristic function of heat

$$\chi\_{\tau}(\nu) = \sqrt{\frac{\det\left(\Lambda\_{\beta z} - i\Lambda\_{\nu z}\right)}{\det\left[\Lambda\_{\beta z} - i\tilde{\Lambda}\_{\nu z}(\tau)\right]}}.\tag{A41}$$

We use the following integral formula:

$$\int d\mathbf{x} e^{-\frac{1}{2}\mathbf{x}^{\mathrm{T}}\mathbf{T}\mathbf{x}} = \sqrt{\frac{(2\pi)^{\dim(\mathbf{T})}}{\det(\mathbf{T})}},\tag{A42}$$

where all the eigenvalues of **T** have positive real parts.

By introducing a diagonal matrix **A** = **<sup>Λ</sup>***βz* − *i***Λ***νz*, we rewrite Equation (A41) as

$$\chi\_{\tau}(\nu) = \sqrt{\frac{1}{\det\left(\mathbf{I} + i\sqrt{\mathbf{A}^{-1}} \left[\boldsymbol{\Lambda}\_{\nu z} - \bar{\boldsymbol{\Lambda}}\_{\nu z}(\tau)\right] \sqrt{\mathbf{A}^{-1}}\right)}}.\tag{A43}$$

Since **Λ ˜** *νz*(*τ*) is a rank-2 matrix, we rewrite it as

$$\tilde{\mathbf{A}}\_{\nu z}(\tau) = 2m\_0 \omega\_0 \tan\left(\frac{\omega\_0 \hbar \nu}{2}\right) \left[ \mathbf{v}\_{q0}(\tau) \mathbf{v}\_{q0}^T(\tau) + \frac{1}{m\_0^2 \omega\_0^2} \mathbf{v}\_{p0}(\tau) \mathbf{v}\_{p0}^T(\tau) \right],\tag{A44}$$

with the vectors

$$\mathbf{v}\_{\theta0}(\mathbf{r}) = \left(\mathfrak{a}\_0, \mathfrak{a}\_1, \dots, \mathfrak{a}\_{N\_r} \frac{\mathfrak{f}\_0}{m\_0}, \frac{\mathfrak{f}\_1}{m\_1}, \dots \frac{\mathfrak{f}\_N}{m\_N}\right)^T,\tag{A45}$$

$$\mathbf{v}\_{p\_0}(\tau) = \left(m\_0 \dot{\mathbf{a}}\_0, m\_0 \dot{\mathbf{a}}\_1, \dots, m\_0 \dot{\mathbf{a}}\_N, \dot{\boldsymbol{\beta}}\_0, \frac{m\_0}{m\_1} \dot{\boldsymbol{\beta}}\_1, \dots, \frac{m\_0}{m\_N} \dot{\boldsymbol{\beta}}\_N\right)^T. \tag{A46}$$

where the evolution time *t* in the terms *αn* and *βn* is set to *τ*. We rewrite the matrix in the determinant (see Equation (A43)) as

$$
\sqrt{\mathbf{A}^{-1}} \left( \Lambda\_{\nu z} - \tilde{\Lambda}\_{\nu z}(\tau) \right) \sqrt{\mathbf{A}^{-1}} = \mathbf{M} \mathbf{M}^{\mathrm{T}},\tag{A47}
$$

with the matrix

$$\mathbf{M}^{\mathrm{T}} = \begin{pmatrix} \sqrt{2m\_{0}\omega\_{0}\tan\frac{\omega\_{0}\hbar\nu}{2}}\mathbf{v}\_{q\_{0}}^{\mathrm{T}}(0) \\ \sqrt{\frac{2}{m\_{0}\omega\_{0}}\tan\frac{\omega\_{0}\hbar\nu}{2}}\mathbf{v}\_{p\_{0}}^{\mathrm{T}}(0) \\ i\sqrt{2m\_{0}\omega\_{0}\tan\frac{\omega\_{0}\hbar\nu}{2}}\mathbf{v}\_{q\_{0}}^{\mathrm{T}}(\tau) \\ i\sqrt{\frac{2}{m\_{0}\omega\_{0}}\tan\frac{\omega\_{0}\hbar\nu}{2}}\mathbf{v}\_{p\_{0}}^{\mathrm{T}}(\tau) \end{pmatrix} \sqrt{\mathbf{A}^{-1}}.\tag{A48}$$

The determinant in Equation (A43) can be simplified to

$$\det\left(\mathbf{I} + i\sqrt{\mathbf{A}^{-1}}\left[\mathbf{A}\_{\nu z} - \tilde{\mathbf{A}}\_{\nu z}(\tau)\right]\sqrt{\mathbf{A}^{-1}}\right) = \det\left(\mathbf{I}\_4 + i\mathbf{M}^T\mathbf{M}\right),\tag{A49}$$

where the right-hand side is the determinant of a 4 × 4 matrix and **I**4 is the 4 × 4 identity matrix. Notice that the initial values of the two vectors are

$$\mathbf{v}\_{\oplus}(0) = \begin{pmatrix} 1, 0, ..., 0, 0, 0, ..., 0 \end{pmatrix}^{T}, \tag{A50}$$

$$\mathbf{v}\_{\mathbb{P}0}(0) = \begin{pmatrix} 0, 0, \ldots, 0, 1, 0, \ldots, 0 \end{pmatrix}^{\mathrm{T}}.\tag{A51}$$

The explicit result of **M**T**M** is obtained as

$$\mathbf{M}^{\mathrm{T}}\mathbf{M} = \boldsymbol{\Xi} \begin{pmatrix} 1 & 0 & i\boldsymbol{\alpha}\_{0} & i\boldsymbol{\alpha}\_{0}/\omega\_{0} \\ 0 & 1 & i\omega\_{0}\boldsymbol{\beta}\_{0} & i\boldsymbol{\beta}\_{0} \\ i\boldsymbol{\alpha}\_{0} & i\omega\_{0}\boldsymbol{\beta}\_{0} & -h\_{11}(\boldsymbol{\tau}) & -h\_{12}(\boldsymbol{\tau}) \\ i\boldsymbol{\alpha}\_{0}/\omega\_{0} & i\boldsymbol{\beta}\_{0} & -h\_{12}(\boldsymbol{\tau}) & -h\_{22}(\boldsymbol{\tau}) \end{pmatrix} \tag{A52}$$

where the elements are functions of the final time *τ* and the functions *h*11(*τ*), *h*12(*τ*) and *h*22(*τ*) are

$$h\_{11}(\tau) = \frac{2m\_0\omega\_0}{\Xi} \tan\left(\frac{\omega\_0\hbar\nu}{2}\right) \mathbf{v}\_{q0}^{\mathrm{T}}(\tau) \mathbf{A}^{-1} \mathbf{v}\_{q0}(\tau),\tag{A53}$$

$$h\_{22}(\tau) = \frac{2}{m\_0 \omega\_0 \Sigma} \tan\left(\frac{\omega\_0 \hbar \nu}{2}\right) \mathbf{v}\_{p\_0}^{\mathrm{T}}(\tau) \mathbf{A}^{-1} \mathbf{v}\_{p\_0}(\tau),\tag{A54}$$

$$\mathcal{H}\_{12}(\tau) = \frac{2}{\Xi} \tan \left( \frac{\omega\_0 \hbar \nu}{2} \right) \mathbf{v}\_{q\_0}^{\mathrm{T}}(\tau) \mathbf{A}^{-1} \mathbf{v}\_{p\_0}(\tau), \tag{A55}$$

with

$$\mathbf{v}\_{q\_0}^T(t)\mathbf{A}^{-1}\mathbf{v}\_{q\_0}(t) = \frac{a\_0^2 + \omega\_0^2\beta\_0^2}{2m\_0\omega\_0\left\{\tanh\left[\frac{(\beta^\* + i\nu)\hbar\omega\_0}{2}\right] - i\tan\left(\frac{\omega\_0\hbar v}{2}\right)\right\}} + \sum\_{n=1}^N \left(\frac{a\_n^2}{\theta\_n} + \frac{1}{\mu\_n}\frac{\beta\_n^2}{m\_n^2}\right),\tag{A56}$$

$$\mathbf{v}\_{p\_0}^T(t)\mathbf{A}^{-1}\mathbf{v}\_{p\_0}(t) = \frac{m\_0\left(\dot{a}\_0^2 + \omega\_0^2\dot{\theta}\_0^2\right)}{2\omega\upsilon\left\{\tanh\left[\frac{(\beta^\*+i\nu)\hbar\omega\_0}{2}\right] - i\tan\left(\frac{\omega\_0\hbar\nu}{2}\right)\right\}} + m\_0^2\sum\_{n=1}^N \left(\frac{\dot{a}\_n^2}{\theta\_n} + \frac{1}{\mu\_n}\frac{\dot{\theta}\_n^2}{m\_n^2}\right),\tag{A57}$$

$$\mathbf{v}\_{q\_0}^T(t)\mathbf{A}^{-1}\mathbf{v}\_{p\_0}(t) = \frac{\frac{d}{dt}\left(a\_0^2 + \omega\_0^2\beta\_0^2\right)}{4\omega\_0\left\{\tanh\left[\frac{(\beta^\*+i\eta)\hbar\omega\_0}{2}\right] - i\tan\left(\frac{\omega\_0\hbar v}{2}\right)\right\}} + \frac{m\_0}{2}\sum\_{n=1}^N \frac{d}{dt}\left(\frac{a\_n^2}{\theta\_n} + \frac{1}{\mu\_n}\frac{\beta\_n^2}{m\_n^2}\right). \tag{A58}$$

The summations are replaced by the integral with the Ohmic spectral density and every element in Equation (A52) is carried out as

$$\alpha\_0(\tau) = \varepsilon^{-\frac{\kappa \tau}{2}} \left[ \cos(\hat{\omega}\_0 \tau) + \frac{\kappa \sin(\hat{\omega}\_0 \tau)}{2\hat{\omega}\_0} \right],\tag{A59}$$

$$\beta\_0(\tau) = \frac{e^{-\frac{k\tau}{2}}\sin(\hat{\omega}v\tau)}{\hat{\omega}\_0},\tag{A60}$$

$$h\_{11}(\tau) = \Theta + e^{-\kappa \tau} \left[ \frac{\omega\_0^2}{\hat{\omega}\_0^2} + \frac{\kappa \sin(2\hat{\omega}\_0 \tau)}{2\hat{\omega}\_0} - \frac{\kappa^2 \cos(2\hat{\omega}\_0 \tau)}{4\hat{\omega}\_0^2} \right] (1 - \Theta), \tag{A61}$$

$$h22(\tau) = \Theta + e^{-\kappa \tau} \left[ \frac{\omega\_0^2}{\hat{\omega}\_0^2} - \frac{\kappa \sin(2\hat{\omega}\_0 \tau)}{2\hat{\omega}\_0} - \frac{\kappa^2 \cos(2\hat{\omega}\_0 \tau)}{4\hat{\omega}\_0^2} \right] (1 - \Theta), \tag{A62}$$

$$h\_{12}(\tau) = \frac{\kappa \omega\_0 e^{-\kappa \tau}}{2\hat{\omega}\_0^2} (\Theta - 1)[1 - \cos(2\hat{\omega}\_0 \tau)].\tag{A63}$$

Then, Equation (17) is obtained by directly calculating the determinant of a 4 × 4 matrix in Equation (A49).

#### **Appendix B. Classical Caldeira–Leggett Model**

We consider the classical Caldeira–Leggett model, where coordinates and momenta commute with each other. To eliminate the initial slippage, the initial state is amended as a coupled state,

$$\rho^{\rm cl}(\mathbf{z};0) = \frac{e^{-\beta' H\_S(0) - \beta[H\_B(0) + H\_{SR}(0)]}}{Z^{\rm cl}(\beta', \beta)},\tag{A64}$$

which represents the probability density in the phase space of the composite system. The classical partition function is obtained by performing the integral in the phase space

$$Z^{\rm cl}(\beta',\beta) = \iint \varepsilon^{-\beta'H\_5(0) - \beta[H\_B(0) + H\_{\rm S\beta}(0)]} dq\_0 dq\_1 \dots dq\_N dp\_0 dp\_1 \dots dp\_N \tag{A65}$$

$$=\frac{2\pi}{\beta'\omega\_0}\prod\_{n=1}^{N}\left(\frac{2\pi}{\beta\omega\_n}\right),\tag{A66}$$

which is independent of the interaction (notice that the partition function of the quantum model relies on the interaction strength [68,87]).

We also define the classical fluctuating heat as the energy difference of the initial and the final system energy. For classical dynamics, the initial and the final states are directly represented by the points in the phase space and the measurements over the system can be applied without disturbing the composite system. Therefore, the characteristic function of heat is

$$\chi^{\rm cl}\_{\tau}(\nu) = \frac{\iint \varepsilon^{j\nu} H\_{\rm S}(\tau) - (\mathcal{J}^{\prime} + i\nu) H\_{\rm S}(0) - \beta[H\_{\rm B}(0) + H\_{\rm S B}(0)]\_{\rm d} dq\_{0} dq\_{1} \dots dq\_{N} dp\_{0} dp\_{1} \dots dp\_{N}}{Z^{\rm cl}(\mathcal{J}^{\prime}, \mathcal{J})},\tag{A67}$$

where the energy of the system *HS*(*t*)=[*p*0(*t*)]2/(<sup>2</sup>*m*0) + *<sup>m</sup>*0*<sup>ω</sup>*20[*q*0(*t*)]2/2 is determined by *p*0(*t*) and *q*0(*t*) associated with the initial point **z**. We choose a new set of initial variables *q*0, q*n* := *qn* − *Cnq*0/(*mnω*2*n*), *p*0 and *pn* in the following calculation.

We rewrite the evolution of the coordinate *q*0(*t*) of the system Equation (A16) as

$$(s^2 + \kappa s + \omega\_0^2)\ddot{q}\_0(s) = \dot{q}\_0(0) + \left(s + \sum\_n \frac{\mathbb{C}\_n^2}{m\_0 m\_n \omega\_n^2} \frac{s}{s^2 + \omega\_n^2} \right) \dot{q}\_0(0) + \sum\_n \frac{\mathbb{C}\_n}{m\_0} \frac{\dot{q}\_n(0) + s\mathbf{q}\_n}{s^2 + \omega\_n^2}.\tag{A68}$$

For the Ohmic spectral density, the summation in the second term is

$$\sum\_{n} \frac{\mathcal{C}\_{n}^{2}}{m\_{0}m\_{n}\omega\_{n}^{2}} \frac{\mathbf{s}}{\mathbf{s}^{2} + \omega\_{n}^{2}} = \boldsymbol{\kappa}.\tag{A69}$$

Thus, Equation (A68) naturally leads to Equation (A19) by substituting *qn*(0) into q*<sup>n</sup>*. The initial slippage is rationally eliminated by choosing a coupled initial state Equation (A64). In reality, the interaction between the system and the heat bath always exists and one cannot prepare the initial state of the composite system without the influence of the interaction. The initial state of the composite system is more likely in the coupled form Equation (A64). The heat bath encodes partial information of the system due to the interaction.

Similar to Equation (A10), the system energy at time *t* can be represented by the dynamical map as

$$H\_{\mathbb{S}}(t) = \frac{1}{2} \mathbf{z}^{\mathrm{T}}(t) \mathbf{A}\_{H\_{\mathbb{S}}} \mathbf{z}(t) \tag{A70}$$

$$\mathbf{z} = \frac{1}{2}\mathbf{z}^{\mathrm{T}}(0)\mathbf{\tilde{A}}\_{H\_{\mathrm{S}}}(t)\mathbf{z}(0),\tag{A71}$$

with the following matrices [88]

$$\mathbf{A}\_{H\_{\mathcal{S}}} = \begin{pmatrix} \frac{m\_0 \omega\_0^2}{2} & & & \\ & \mathbf{O} & & \\ & & \frac{1}{m\_0} & \\ & & & \mathbf{O} \end{pmatrix} \,\tag{A72}$$

and

$$\tilde{\mathbf{A}}\_{H\_{\rm S}}(t) = \exp(\mathbf{L}^{\rm T}t)\mathbf{A}\_{H\_{\rm S}}\exp(\mathbf{L}t). \tag{A73}$$

The initial vector is now amended to

$$\mathbf{z}(0) = \left(\mathfrak{q}\_0(0), \mathfrak{q}\_1, \dots, \mathfrak{q}\_{\mathfrak{u}}, p\_0(0), \dots, p\_N(0)\right)^T. \tag{A74}$$

The initial Hamiltonians *HS*(0) and *HB*(0) + *HSB*(0) are

$$H\_S(0) = \frac{1}{2} \mathbf{z}^T(0) \mathbf{A}\_{H\_S} \mathbf{z}(0),\tag{A75}$$

$$\begin{split} H\_B(0) + H\_{SB}(0) &= \sum\_{n=1}^{N} \left( \frac{1}{2} \frac{p\_n^2}{m\_n} + \frac{1}{2} m\_n \omega\_n^2 \mathbf{q}\_n^2 \right) \\ &= \frac{1}{2} \mathbf{z}^{\mathrm{T}}(0) \mathbf{A}\_{H\_B} \mathbf{z}(0), \end{split} \tag{A76}$$

with the matrix

$$\mathbf{A}\_{H\_B} = \text{diag}(0, m\_1 \omega\_{1'}^2, \dots, m\_N \omega\_{N'}^2, 0, \frac{1}{m\_1}, \dots, \frac{1}{m\_N}). \tag{A77}$$

According to the integral Formula (A42), we carry out the characteristic function Equation (A67) into

$$\chi\_{\tau}^{\rm cl}(\nu) = \sqrt{\frac{\det\left[\beta^{\prime}\mathbf{A}\_{H\_{\rm S}} + \beta\mathbf{A}\_{H\_{\rm B}}\right]}{\det\left[\beta^{\prime}\mathbf{A}\_{H\_{\rm S}} + \beta\mathbf{A}\_{H\_{\rm B}} - i\nu\left(\tilde{\mathbf{A}}\_{H\_{\rm S}}(\tau) - \mathbf{A}\_{H\_{\rm S}}\right)\right]}}.\tag{A78}$$

For the classical limit ¯*h* → 0, we can verify

$$\lim\_{\hbar \to 0} \frac{\mathbf{A}\_{\beta z}}{\hbar} = (\beta' + i\nu)\mathbf{A}\_{H\_{\mathbb{S}}} + \beta \mathbf{A}\_{H\_{\mathbb{B}'}} \tag{A79}$$

$$\lim\_{\hbar \to 0} \frac{\mathbf{A}\_{\nu z}}{\hbar} = \nu \mathbf{A}\_{H\_{\rm S}} \tag{A80}$$

$$\lim\_{\hbar \to 0} \frac{\ddot{\Lambda}\_{\nu z}(t)}{\hbar} = \nu \ddot{\Lambda}\_{H\_{\S}}(t),\tag{A81}$$

and obtain

$$\lim\_{\hbar \to 0} \chi\_{\mathbb{T}}(\nu) = \chi\_{\mathbb{T}}^{\mathrm{cl}}(\nu),\tag{A82}$$

with *χτ*(*ν*) given in Equation (A41). The final result Equation (A78) is the same as Equation (25). Hence, we use the same notation.

#### **Appendix C. The Characteristic Function of Heat for the Classical Brownian Motion**

We derive the characteristic function of heat for the classical Brownian motion. For an underdamped Brownian particle moving in a potential *<sup>V</sup>*(*x*), the stochastic dynamics is described by the complete Langevin equation

$$
\ddot{\mathbf{x}} + \kappa \dot{\mathbf{x}} + \frac{1}{m} \frac{\partial V}{\partial \mathbf{x}} = \frac{1}{m} F\_{\text{fluc}}(t). \tag{A83}
$$

The fluctuating force *<sup>F</sup>*fluc(*t*) is a Gaussian white noise satisfying the fluctuation–dissipation relation

$$
\langle F\_{\rm fluc}(t) F\_{\rm fluc}(t') \rangle = 2m\kappa k\_B T \delta(t - t'). \tag{A84}
$$

The evolution of the system state is characterized by the probability density function *ρ*(*<sup>x</sup>*, *p*; *t*) in the phase space. The stochastic dynamics is then described by the Kramers equation [86]

$$\frac{\partial \rho}{\partial t} = \mathcal{L}^{\rho}[\rho],\tag{A85}$$

with the Liouville operator

$$\mathcal{L}^{\rho}[\rho] = -\frac{\partial}{\partial \mathbf{x}}(\frac{p}{m}\rho) + \frac{\partial}{\partial p}\left[\kappa p \rho + \frac{\partial V(\mathbf{x})}{\partial \mathbf{x}}\rho + \frac{\kappa m}{\beta}\frac{\partial \rho}{\partial p}\right]. \tag{A86}$$

Similarly, in the phase space, we calculate the characteristic function of heat for the classical Brownian motion

$$\chi^{\rm cl}\_{\tau}(\nu) = \iint d\mathbf{x}d\boldsymbol{p} \varepsilon^{\rm iv} \left[ \frac{p^2}{2\overline{u}} + V(\boldsymbol{x}) \right] \boldsymbol{\eta} (\mathbf{x}, \boldsymbol{p}; \tau), \tag{A87}$$

where a probability-density-like function *η*(*<sup>x</sup>*, *p*; *t*) also satisfies the dynamic Equation (A85) with the initial condition

$$\eta(\mathbf{x}, p; 0) = e^{-i\nu \left[\frac{p^2}{2m} + V(\mathbf{x})\right]} \rho(\mathbf{x}, p; 0). \tag{A88}$$

We consider the system potential as a harmonic potential *<sup>V</sup>*(*x*) = *mω*20 *x*2/2 and the initial system state as an equilibrium state

$$\rho(\mathbf{x}, p; \mathbf{0}) = \frac{1}{Z\_S^{\rm cl}(\boldsymbol{\beta}')} e^{-\boldsymbol{\beta}' \left(\frac{\boldsymbol{p}^2}{2\pi} + \frac{1}{2}m\omega\_0^2 \mathbf{x}^2\right)\_\star} \tag{A89}$$

with the inverse temperature *β* and the classical partition function *Z*cl *S* (*β*) = <sup>2</sup>*π*/(*βω*0). Under such conditions, the probability-density-like function *η*(*<sup>x</sup>*, *p*; *t*) is always in a quadratic form, assumed as

$$\eta(\mathbf{x}, p; t) = \frac{1}{Z\_{\mathcal{S}}^{\rm cl}(\beta')} \mathbf{c}^{-\left[a(t)\frac{p^2}{2m} + b(t)\frac{1}{2}m\omega\_0^2 \mathbf{x}^2 + c(t)\omega\_0 \mathbf{x} p + \Lambda(t)\right]} \,. \tag{A90}$$

The Kramers Equation (A85) for *η*(*<sup>x</sup>*, *p*; *t*) leads to the following ordinary differential equations:

$$
\Lambda = -\kappa \left( 1 - \frac{a}{\beta} \right),
\tag{A91}
$$

$$
\dot{a} = 2\kappa a \left( 1 - \frac{a}{\beta} \right) - 2\omega\_0 c,\tag{A92}
$$

$$
\dot{b} = 2c \left( \omega\_0 - \frac{\kappa}{\beta} c \right),
\tag{A93}
$$

$$
\dot{\varepsilon} = \omega \upsilon (a - b) + \kappa c - 2 \frac{\kappa}{\beta} ac\_\prime \tag{A94}
$$

with the initial conditions *a*(0) = *b*(0) = *β* + *i<sup>ν</sup>*, *c*(0) = 0 and Λ(0) = 0. According to the conservation of the probability *η*(*<sup>x</sup>*, *p*; *t*)*dxdp* = const, the coefficient Λ(*t*) is obtained as

$$e^{-\Lambda(t)} = \frac{\sqrt{a(t)b(t) - c(t)^2}}{\beta' + i\nu}. \tag{A95}$$

Substituting Equation (A90) into Equation (A87), we obtain the characteristic function for the classical Brownian motion as

$$\chi^{\rm cl}\_{\tau}(\nu) = \frac{\beta'}{\beta' + i\nu} \sqrt{\frac{a(\tau)b(\tau) - c(\tau)^2}{[a(\tau) - i\nu][b(\tau) - i\nu] - c(\tau)^2}}. \tag{A96}$$

To solve the nonlinear differential Equations (A92)–(A94), we introduce a new set of variables,

$$A = \frac{a}{ab - c^2},\tag{A97}$$

$$B = \frac{b}{ab - c^2},\tag{A98}$$

$$\mathcal{C} = \frac{c}{ab - c^2},\tag{A99}$$

and obtain the linear differential equations

$$\frac{dA}{dt} = -2\omega\rho\mathcal{C},\tag{A100}$$

$$\frac{dB}{dt} = 2\omega\_0 \mathbf{C} - 2\mathbf{x}B + 2\frac{\kappa}{\beta} \mathbf{y} \tag{A101}$$

$$\frac{d\mathbb{C}}{dt} = \omega\_0 (A - B) - \kappa \mathbb{C},\tag{A102}$$

with the initial conditions *A*(0) = *B*(0) = 1/(*β* + *<sup>i</sup>ν*) and *C*(0) = 0. The characteristic function Equation (A96) becomes

$$\chi\_{\rm tr}^{\rm cl}(\nu) = \frac{\beta'}{\beta' + i\nu} \sqrt{\frac{1}{1 - i\nu(A + B) - \nu^2(AB - C^2)}}.\tag{A103}$$

The solutions to Equations (A100)–(A102) are

$$A(t) = -\frac{e^{-\kappa t} \left[\kappa^2 \cos(2\hat{\omega}\_0 t) - 2\hat{\omega}\_0 \kappa \sin(2\hat{\omega}\_0 t) - 4\omega\_0^2\right]}{4\beta\hat{\omega}\_0^2} \frac{\not\!\!/ - \not\!\!/ - i\nu}{\not\!\!/ ' + i\nu} + \frac{1}{\beta},\tag{A104}$$

$$B(t) = -\frac{e^{-\kappa t} \left[\kappa^2 \cos(2\omega\_0 t) + 2\omega\_0 \kappa \sin(2\omega\_0 t) - 4\omega\_0^2\right]}{4\beta\omega\_0^2} \frac{\not\!\!/ - \not\!\!/ - i\nu}{\not\!\!/ ' + i\nu} + \frac{1}{\beta},\tag{A105}$$

$$\mathcal{L}(t) = -\frac{2e^{-\mathbf{x}t}\kappa\omega\_0[\cos(2\hat{\omega}\_0 t) - 1]}{4\beta\hat{\omega}\_0^2} \frac{\beta - \beta' - i\nu}{\beta' + i\nu},\tag{A106}$$

with *ω*<sup>ˆ</sup> 0 = *ω*<sup>20</sup> − *κ*2/4. Plugging the solutions into Equation (A103), we immediately obtain Equation (25). We remark that the heat distribution of the classical Brownian motion was obtained by the path-integral method in Ref. [58], but they only consider the initial temperature of the system to be the same as that of the bath.

#### *Appendix C.1. The Long-Time Limit*

After sufficiently long relaxation time, the solutions *<sup>a</sup>*(*t*), *b*(*t*) and *c*(*t*) to Equations (A92)–(A94) eventually approach *a*(∞) = *b*(∞) = *β* and *c*(∞) = 0. The longtime limit of Equation (A96) reproduces Equation (33).

#### *Appendix C.2. The Underdamped Limit*

In the underdamped limit *<sup>κ</sup>*/*<sup>ω</sup>*0 → 0, the differential Equations (A92)–(A94) are reduced to

$$
\dot{a} = \kappa a \left( 1 - \frac{a}{\beta} \right),
\tag{A107}
$$

with *b* = *a* and *c* = 0. The solution is

$$a(t) = \frac{\beta(\beta' + i\nu)}{\beta' + i\nu + (\beta - \beta' - i\nu)e^{-\kappa t}},\tag{A108}$$

and Equation (A96) becomes

$$\chi\_{\tau}^{\text{w,cl}}(\nu) = \frac{\beta \beta'}{(\beta - i\nu)(\beta' + i\nu)(1 - e^{-\kappa \tau}) + \beta \beta' e^{-\kappa \tau}}. \tag{A109}$$

It can be checked that Equation (37) reproduces Equation (A109) in the classical limit *h*¯ → 0.

*Appendix C.3. The Overdamped Limit*

In the overdamped limit *<sup>κ</sup>*/*<sup>ω</sup>*0 → <sup>∞</sup>, the differential Equations (A92)–(A94) are reduced to

$$
\dot{a} = 2\kappa a \left( 1 - \frac{a}{\beta} \right),
\tag{A110}
$$

$$
\dot{\theta} = 2\omega\wp c \left( 1 - \frac{\kappa}{\beta\omega\_0}c \right),
\tag{A111}
$$

$$0 = (\kappa c + \omega\_0 a) - 2\frac{\kappa}{\beta}ac - b\omega\_0. \tag{A112}$$

Eliminating *c* in Equation (A112), Equation (A111) becomes

$$\dot{\beta} = \frac{2\omega\_0^2}{\kappa} \frac{(\beta - a - b)(b - a)}{(\beta - 2a)\left(1 - 2\frac{a}{\beta}\right)}.\tag{A113}$$

Notice that, in the overdamped limit, the relaxation timescales of the momentum and the coordinate are separated. We can substitute *a* = *β* in Equation (A113) and obtain

$$
\beta = \frac{2\omega\_0^2}{\kappa} b \left( 1 - \frac{b}{\beta} \right). \tag{A114}
$$

With the initial condition *a*(0) = *β* + *iν* and *b*(0) = *β* + *i<sup>ν</sup>*, the solutions are

$$a(t) = \frac{(\beta' + i\nu)\beta}{\beta' + i\nu + (\beta - \beta' - i\nu)e^{-2\pi t}}\tag{A115}$$

$$b(t) = \frac{(\beta' + i\nu)\beta}{\beta' + i\nu + (\beta - \beta' - i\nu)e^{-\frac{2\omega\_0^2}{\kappa}t}}.\tag{A116}$$

We substitute Equations (A115), (A116) and *c*(*t*) ≈ 0 into Equation (A96) and obtain

$$\begin{split} \chi^{\text{scl}}\_{\tau}(\nu) &= \frac{\beta \beta'}{\sqrt{(\beta - i\nu)(\beta' + i\nu)(1 - e^{-2\kappa\tau}) + \beta\beta'e^{-2\kappa\tau}}} \\ &\times \frac{1}{\sqrt{(\beta - i\nu)(\beta' + i\nu)\left(1 - e^{-\frac{2\omega\_0^2}{\kappa}\tau}\right) + \beta\beta'e^{-\frac{2\omega\_0^2}{\kappa}\tau}}}. \end{split} \tag{A117}$$

It can be checked that Equation (38) reproduces Equation (A117) in the classical limit *h*¯ → 0.

#### **Appendix D. The Characteristic Function of Heat for Complete Thermalization**

We derive the characteristic function of heat for a complete thermalization process, Equation (34), in the main content. For a complete thermalization process (typically with infinite relaxation time), the information of the initial state is completely forgotten and the final state is always an equilibrium state at the inverse temperature *β* of the heat bath,

$$
\gamma\_{\text{th}\downarrow^{\nu\_I}} = p\_{l'}^{\text{eq}},\tag{A118}
$$

regardless of the initial state *l*. Therefore, the characteristic function of heat for complete thermalization is *<sup>χ</sup>*th(*ν*) = ∑*l*,*<sup>l</sup>* exp[*<sup>i</sup>ν*(*ESl* − *ESl* )]*pl <sup>p</sup>*eq*l* . We immediately obtain Equation (34) by plugging into the initial distribution *pl* = exp(−*βESl* )/*ZS*(*β*) and the final distribution *<sup>p</sup>*eq*l* = exp(−*βESl*)/*ZS*(*β*), where *ZS*(*β*) = ∑*l* exp(−*βESl* ) is the partition function of the system at the inverse temperature *β*.

#### **References and Notes**



#### *Article* **Does Decoherence Select the Pointer Basis of a Quantum Meter?**

**Abraham G. Kofman \* and Gershon Kurizki \***

> Department of Chemical and Biological Physics, Weizmann Institute of Science, Rehovot 761001, Israel

**\*** Correspondence: kofmana@gmail.com or abraham.kofman@weizmann.ac.il (A.G.K.);

gershon.kurizki@weizmann.ac.il (G.K.)

**Abstract:** The consensus regarding quantum measurements rests on two statements: (i) von Neumann's standard quantum measurement theory leaves undetermined the basis in which observables are measured, and (ii) the environmental decoherence of the measuring device (the "meter") unambiguously determines the measuring ("pointer") basis. The latter statement means that the environment *monitors* (measures) *selected* observables of the meter and (indirectly) of the system. Equivalently, a measured quantum state must end up in one of the "pointer states" that persist in the presence of the environment. We find that, unless we restrict ourselves to projective measurements, decoherence does not necessarily determine the pointer basis of the meter. Namely, generalized measurements commonly allow the observer to choose from a multitude of alternative pointer bases that provide the same information on the observables, regardless of decoherence. By contrast, the measured observable does not depend on the pointer basis, whether in the presence or in the absence of decoherence. These results gran<sup>t</sup> further support to our notion of Quantum Lamarckism, whereby the observer's choices play an indispensable role in quantum mechanics.

**Keywords:** quantum measurements; decoherence; pointer states; Quantum Lamarckism; the observer in quantum mechanics
