Asymmetries of thermal processes in open quantum systems

Abstract

An intriguing phenomenon in non-equilibrium quantum thermodynamics is the asymmetry of thermal processes. Relaxation to thermal equilibrium is the most important dissipative process, being a key concept for the design of heat engines and refrigerators, contributing to the study of foundational questions of thermodynamics, and being relevant for quantum computing through the process of algorithmic cooling. Despite the importance of this kind of processes, their dynamics are far from being understood. We show that the free relaxation to thermal equilibrium follows intrinsically different paths depending on whether the temperature of the system increases (heating up) or decreases (cooling down), being faster in the first case. Our theory is exemplified using the recently developed thermal kinematics based on information geometry theory, utilizing three prototypical examples: a quantum two-level system, the quantum harmonic oscillator, and a trapped quantum Brownian particle, including both analytic results and numerical simulations. For this, we have extended the thermal kinematic approach to open quantum systems. Additionally, we offer a simple theoretical explanation in the case of a two level system and a more general picture for the other two systems based on the spectral decomposition of the Liouvillian and the spectral gap of reciprocal processes.

I Introduction

When a system is pushed far from equilibrium, its evolution may follow anomalous paths. A series of seminal works done during the past century [1, 2, 3, 4, 5, 6, 7] has provided essential advances in studying transitory phenomena in the linear regime associated with fluctuations, except for some particular cases [8, 9] where predictions can extend beyond equilibrium. Despite this progress, we still lack a general theory beyond linear response and fluctuation theorems to decipher the dynamics and behavior of transient regimes of a freely evolving system between two desired states [10, 11]. This problem is of particular interest for quantum thermodynamic processes [12], finite-time quantum heat engines [13, 14, 15, 16] and establishing speed limit bounds [17, 18, 19, 20]. Recent progress in unraveling anomalous shortcuts during relaxation processes in out-of-equilibrium systems points in this direction [21].

A remarkable example of a possible counter-intuitive behavior of a system is the Mpemba-like effect (ME) [22, 23, 24]. Namely, put two identical systems at different initial temperatures in contact with a reservoir at a hotter or colder temperature than those of the two systems. The ME occurs when the initially hotter/colder system cools down/heats up faster than the system that was initially closer to the final temperature. In the case of cooling, the effect is called normal ME, and for heating, it is called inverse ME [25, 26]. In Markovian systems, the ME can be well understood using a spectral decomposition of the decay modes, diminishing (weak ME), or canceling slow-decaying modes (strong ME) to enhance the fast ones, making it possible to control the speed of the relaxation. In this way, up to an exponential acceleration is achievable [27]. This phenomenon has been realized both in classical [25, 27, 28, 29, 30, 31, 32, 33] and open quantum systems [34, 35, 36, 37, 38]. Additionally, a generalization of the ME to quantum entangled configurations has been very recently proposed [39, 40, 41, 42]. Note that a strong relation exists between exceptional points and speed up relaxation in open quantum systems [36, 43].

Alternatively, when spectral methods are not applicable, other strategies can be used to understand anomalous evolution using macroscopic observables depending on the system of interest. The origin of anomalous relaxation is associated with energy non-equipartition in water and granular gases composed by rough hard spheres [44, 45], a particular condition in kurtosis also in the former with smooth hard spheres [26, 46], and correlation length in spin glasses [47]. Furthermore, the strategy of employing several sudden changes in temperature has been probed useful for shortening relaxation times, such as preheating protocols [48]. This approach takes advantage of the slow growth of magnetic domains near phase transitions in systems where time-scale separation is not possible [49], or through different control techniques [50, 51, 52].

Figure 1: (a) Asymmetric cooling and heating relaxation to an equidistant stationary state at temperature $T_{W}$. It takes longer when the system is initially hot (thermalized with a bath at $T_{H}$ which is decoupled at $t=0$) than when it is cold (at $T_{C}$), with $T_{C}<T_{W}<T_{H}$. (b) Asymmetric cooling and heating evolution between two states at temperatures $T_{C}$ and $T_{H}$, with $T_{C}<T_{H}$. The evolution from hot to temperature (cooling) is slower than from cold to hot (heating).

A fundamental question, illustrated in Fig. 1, is whether free cooling and heating processes after a sudden change of the environment temperature are identical or follow intrinsically different paths, see Fig. 1(a). In classical systems heating and cooling can show an asymmetry that has been verified both theoretically and experimentally far from equilibrium  [53, 54]. An even more emphatic result is that the asymmetry is revealed when relaxation processes occur between two fixed temperatures [54], see Fig. 1(b). This has been successfully explained mathematically by using the so-called thermal kinematics [54], based on information geometrical arguments [55, 56]. In this paper we focus on that question, that is, unraveling the mechanism of the heating and cooling processes in the realm of open quantum systems. In order to do this we develop an extension of the thermal kinematics theory to the thermodynamics of open quantum systems. We analyze whether a relaxation process far from equilibrium, say from an initially hot to a colder thermal state, is equally fast than its reverse, from the colder to the hotter, and relate it to the properties of the spectral gap [57, 58, 59, 60, 61, 62]. To showcase this, we use simple models based on a thermal qubit, a quantum harmonic oscillator, and a quantum Brownian particle.

The heat properties of such simple quantum systems have recently become accessible experimentally. Solid state realizations of qubits coupled to fermionic or bosonic reservoirs allow to control the spectral properties, couplings and temperatures externally [63, 64, 65]. This is the case of quantum dot systems [66, 67], which can selectively be (un)connected to different reservoirs with gate voltages [68] and whose distribution can be measured via charge detectors [69, 70, 71, 72], or of superconducting circuits coupled to resistors acting as thermal baths via tunable resonators [73, 74, 75, 76, 77, 78]. Furthermore, the qubit state can be monitored [79, 80, 81, 82, 83]. Improvements in high frequency thermometry even allows to detect single temperature fluctuations [84]. These ingredients make the detection of relaxation paths in quantum information systems possible.

The recent measurement of asymmetric relaxation of a classical particle in a harmonic trap [54] motivates us to treat this problem from a quantum perspective. To do so, we investigate the thermalization of a quantum Brownian particle, a model that has successfully been applied to describe a plethora of quantum effects, such as quantum dissipation [85, 86], harmonic oscillators [87], macroscopic quantum tunneling [88, 89, 90], metastable states [91], single-electron transistors [92], the spin-boson problem [93], or impurity dynamics in Luttinger liquids [94] and ultracold atomic gases [95]. We hence emphasize that understanding the relaxation processes is of importance for quantum thermodynamics and for the physics of driven nanoscale devices [12, 11].

This manuscript is organized as follows. In Secs. II and III, we present the theoretical framework based on the master equation for open quantum systems and the measures of thermodynamic distances. In Sec. IV, the definitions of the different protocols are provided. We then apply these methods to three different systems. In Sec. V we consider a the simplest case of two level system coupled to a thermal bath, which can be solved analytically. Then Sec. VI considers more complex systems, namely the harmonic oscillator and the quantum Brownian particle. Additionally, this section includes numerical simulations for the non-analytically solvable systems. Section VII provides a theoretical justification for all the phenomena based on the spectrum of the Liouvillian and the influence of the considered initial state on the evolution. Finally, the conclusions are drawn in Sec. VIII.

II Theoretical framework: Markovian Open Quantum Systems

The state of the quantum system, weakly coupled to the environment, is described by its reduced density matrix $\rho(t)$, whose evolution is governed by the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) quantum master equation [96, 97, 98, 99, 100, 101, 102, 103, 104]

$$\dot{\rho}(t)=\mathcal{L}[\rho(t)],$$ (1)

where $\mathcal{L}$ is the Liouvillian superoperator

$$\mathcal{L}[\rho(t)]=-\dfrac{i}{\hbar}[H,\rho(t)]+\sum_{i=1}^{N}\left(L_{i}\rho(t)L^{\dagger}_{i}-\frac{1}{2}\left\{L^{\dagger}_{i}L_{i},\rho(t)\right\}\right),$$ (2)

being $H$ the Hamiltonian of the system, describing its coherent dynamics. The $N$ jump operators $L_{i}$ describe the dissipative effects due to the presence of an environment. The Liouvillian superoperator $\mathcal{L}$ preserves the trace, i.e. $\Tr\left(\mathcal{L}[\rho(t)]\right)=0$, hermiticity, i.e. $\left(\mathcal{L}[\rho(t)]\right)^{\dagger}=\mathcal{L}[\rho^{\dagger}(t)]$, $\forall\rho(t)$, and complete positivity.

The general solution to Eq. (1) can be directly obtained as $\rho(t)=e^{t\mathcal{L}}[\rho(0)]$, where the superoperator $e^{t\mathcal{L}}$ is defined by its power expansion. Assuming the generator to be diagonalizable, one can find the right eigenmatrices, $\Lambda^{r}_{k}$, such that

$$\mathcal{L}[\Lambda^{r}_{k}]=\lambda_{k}\Lambda^{r}_{k}.$$ (3)

The complex numbers $\lambda_{k}$ are the eigenvalues of the Liouvillian. Note that, due to the hermiticity-preservation of $\mathcal{L}$, if $\lambda_{k}$ is a complex eigenvalue, then $\lambda_{k}^{*}$ must also be an eigenvalue. For the same reason, one can also show that if $\lambda_{k}$ is real, then $\Lambda^{r}_{k}$ can be chosen to be Hermitian. Associated with the map defined in Eq. (2), there is a dual map, also called the adjoint Lindblad map, which implements the evolution of observables:

$$\mathcal{L}^{\dagger}[O]=\dfrac{i}{\hbar}[H,O]+\sum_{i=1}^{N}\left(L^{\dagger}_{i}OL_{i}-\frac{1}{2}\left\{O,L^{\dagger}_{i}L_{i}\right\}\right).$$ (4)

This dual map, $\mathcal{L}^{\dagger}$, is diagonalized by the left eigenmatrices $\Lambda^{\ell}_{k}$,

$$\mathcal{L}^{\dagger}[\Lambda^{\ell}_{k}]=\lambda_{k}\Lambda^{\ell}_{k}.$$ (5)

The matrices $\Lambda^{\ell}_{k}$ are in principle different from the matrices $\Lambda^{r}_{k}$ in Eq. (3). However, $\Lambda^{\ell}_{k}$ and $\Lambda^{r}_{k}$ still form a bi-orthogonal basis for the space of matrices and can always be defined fulfilling the property $\Tr\left(\Lambda^{\ell}_{k}\Lambda^{r}_{h}\right)=\delta_{kh}$.

Since the dynamics generated by $\mathcal{L}$ is completely positive, the eigenvalues of the Liouvillian superoperator all have a non-positive real part, ${\rm Re}\left(\lambda_{k}\right)\leq 0$. Furthermore, for bounded systems, Evan’s Theorem [105] enforces that at least one eigenvalue is zero, $\lambda_{1}=0$, and this is also the case for many unbounded systems. Assuming that the null eigenvalue is non-degenerate, the asymptotic stationary state of the open quantum system is directly related to its associated eigenmatrix [106, 107],

$$\rho_{\rm ss}=\lim_{t\to\infty}\rho(t)=\Lambda_{1}^{r}.$$ (6)

Integrating Eq. (2), the spectral decomposition of $\mathcal{L}$ allows us to write the dynamics of any initial density matrix as

$$\rho(t)=e^{t\mathcal{L}}\left[\rho_{0}\right]=\Lambda^{r}_{1}+\sum_{k=2}^{d^{2}}e^{t\lambda_{k}}{\rm Tr}\left(\Lambda^{\ell}_{k}\rho_{0}\right)\Lambda^{r}_{k},$$ (7)

where $d$ is the dimension of the Hilbert space of the system. This decomposition shows that the matrices $\Lambda_{k}^{r}$ are nothing but the excitation modes of the system, each one characterized by a decay rate $|{\rm Re}(\lambda_{k})|$. For long times, the relevant terms are those related to the $\lambda_{k}$ with the smallest real part in modulus and finite overlap with the initial state. To study the time-evolution of our systems we order the eigenvalues $\lambda_{k}$ in such a way that $|{\rm Re}\left(\lambda_{2}\right)|\leq|{\rm Re}\left(\lambda_{3}\right)|\leq\ldots\leq|{\rm Re}\left(\lambda_{m}\right)|$. The overlap between the $i-$th eigenmatrix and the initial state, $\rho_{0}$, is determined by

$$\xi_{i}=\Tr\left(\Lambda_{i}^{\ell}\rho_{0}\right),$$ (8)

Note that this term $\xi_{i}$ is the same as the one appearing in the sum presented in Eq. (7). This term will provide us with the influence of the Lindbladian, which fixes the temporal evolution, onto the initial state.

III Quantum Thermal kinematics: Measures of distance and speed

The concept of thermal kinematics, well-established for classical systems recently in Ref. [54], combines arguments from stochastic thermodynamics with information geometry to analyze the thermodynamical processes [56]. For classical systems, it is possible to define a statistical distance [54], related to the classical Fisher information which quantifies the temporal variation of local flows. Therefore, for two time-varying infinitesimal processes, the line element can be defined from the Kullback-Leibler divergence (KLD) of two probability distributions, defined as

$$D_{\text{cl}}\left[P_{\text{cl}}(x,t+dt),P_{\text{cl}}(x,t)\right]=I_{\text{cl}}(t)dt^{2}+\mathcal{O}(dt^{4}),$$ (9)

where $I_{\text{cl}}(t)$ is the classical Fisher information (FI), and allows us to define a proper statistical distance between two states (see App. A, Eq. (57)). Note that we denote all classical quantities and variables with the subscript $\mathrm{cl}$. The line element is then defined from Eq. (9) as

$$dl_{\mathrm{cl}}:=\sqrt{I_{\text{cl}}(t)}dt.$$ (10)

where $\sqrt{I_{\text{cl}}(t)}$ can be identified as the statistical velocity at a given time $t$, namely

$$v_{\text{cl}}(t):=\sqrt{I_{\text{cl}}(t)}.$$ (11)

To study thermal kinematics in the quantum regime, we may use two different measures. The first one will be the fidelity between two states (analog to the KLD in the classical case), defined as

$$F(\rho_{1},\rho_{2}):=\Tr\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{1}}}.$$ (12)

It measures how close two quantum states are in terms of their density matrix. It is symmetric and invariant under unitary operations. Despite it does not define a distance [108], the fidelity allows us to define a proper metric, the so-called Bures distance

$$\left[D_{B}(\rho,\sigma)\right]^{2}:=2(1-F(\rho,\sigma)).$$ (13)

Analogous to the classical case, in our context of thermal relaxation, an infinitesimal statistical line element may be defined as follows [109]

$$\left[D_{B}\left(\rho(t),\rho(t+dt)\right)\right]^{2}=\frac{1}{4}\mathcal{I}_{Q}[\rho(t)]dt^{2}+\mathcal{O}(dt^{4}),$$ (14)

being $\mathcal{I}_{Q}$ the quantum Fisher information (QFI), with respect to the parameter time, defined as

$$\mathcal{I}_{Q}[\rho(t)]:=\text{Tr}\left[L_{t}^{2}\rho(t)\right],$$ (15)

where $L_{t}$ is the logarithmic time-derivative operator defined by $\dot{\rho}(t):=\left(L_{t}\rho(t)+\rho(t)L_{t}\right)/2$, see App. A. From Eq. (14), we can directly define the line element as

$$dl:=\sqrt{\frac{1}{4}\mathcal{I}_{Q}[\rho(t)]}dt,$$ (16)

and thus

$$v(t):=\sqrt{\frac{1}{4}\mathcal{I}_{t}[\rho(t)]}$$ (17)

represents the quantum instantaneous statistical velocity of the system in the quantum case. The statistical length of a path taken between time $t_{i}$ and $t_{f}$ is computed as

$$\ell(t_{i},t_{f})=\int_{t_{i}}^{t_{f}}\sqrt{\frac{1}{4}\mathcal{I}_{t}\left[\rho(t)\right]}dt.$$ (18)

As reaching the steady state during a dissipative process takes infinite time, to establish a kinematic basis for quantifying thermal relaxation kinematics, we define the quantum degree of completion as

$$\varphi(s):=\dfrac{\ell(t_{i},t_{s})}{\ell(t_{i},t_{f})},$$ (19)

being a monotonically increasing function bounded between 0 and 1.

IV Heating and cooling protocols

To puzzle out the properties of cooling and heating far from equilibrium in quantum systems subject to instantaneous quenches, we define two possible experiments.

IV.1 Protocol between three temperatures

The first feasible protocol is to compare the free evolution with respect to an intermediate temperature. Hence, we define three temperatures $T_{C}<T_{W}<T_{H}$, the subscripts corresponding to cold ($C$), warm ($W$) and hot ($H$) respectively, as displayed in Fig. 1(a). Associated to this temperatures there are three Gibbs states, $\rho^{\text{th}}_{\beta_{i}}=\exp[-\beta_{i}H]/\mathcal{Z}$, with $H$ being the Hamiltonian of the system, $\beta_{i}=1/k_{B}T_{i}$ the inverse temperature, and $\mathcal{Z}=\Tr\left\{\exp[-\beta_{i}H]\right\}$ the partition function for $i\in\{C,W,H\}$.

In this protocol, we can analyze the behavior of our system by the use of the fidelity, Eq. (12). As both trajectories, cooling and heating up, have the same target point, the warm thermal density matrix, we can use the fidelity between our temporal state and the target one as a measure of distance. To fix the initial conditions, we consider thermal states with equal fidelity values with respect to $T_{W}$ for both $T_{C}$ and $T_{H}$, meaning that

$$F(\rho^{\text{th}}_{\beta_{C}},\rho^{\text{th}}_{\beta_{W}})=F(\rho^{\text{th}}_{\beta_{H}},\rho^{\text{th}}_{\beta_{W}}).$$ (20)

The process starts with the hot and cold temperature and introduce a sudden quench to the warm temperature, monitoring the fidelity evolution. The three temperatures protocol allows us to distinguish which process is faster heating or cooling.

We first focus on what we call forward protocol where the relaxation occurs to the warm temperature, $T_{W}$ between the state at hot temperature $T_{H}$, and cold temperature, $T_{C}$. What may be even more surprising is that heating also turns out to be faster along the reversed, the so called backward protocol. That is, we prepare the system to be in equilibrium at the warm temperature $T_{W}$ and track back the relaxations at $T_{C}$ and $T_{H}$, respectively. This observation is remarkable as it shows that heating is inherently faster than cooling at TE conditions.

IV.2 Protocol between two temperatures

We can also proceed using a more simple protocol, namely, cooling and heating between two temperatures $T_{C}<T_{H}$ respectively, as displayed in Fig. 1(b). Measuring with this second protocol the asymmetry between cooling and heating far from equilibrium appears more resounding and stronger. But in this protocol the absence of a reference density matrix prevents us to use the fidelity as a distance measure. We need to use a true metric distance, namely, the quantum Fisher information, Eq. (15), and the so called thermal kinematics [54].

In this scenario, starting from one of the temperatures, after a sudden quench, we let the system evolve freely to the other one. This phenomenon allows us to observe heating, i.e relaxation at $T_{H}$ in a temperature quench from an equilibrium prepared at $T_{C}$; and the reverse cooling, i.e. relaxation at $T_{C}$ in a temperature quench from the equilibrium at $T_{H}$. In order to compare the two processes in a proper way we will use the quantum degree of completion given by Eq. (19) and the quantum instantaneous statistical velocity, Eq. (17).

In the following, we test these protocols in three different quantum systems of increasing complexity. The first one is the simplest case, as it only consists of a two-level system coupled to a thermal bath at a given temperature. In this case, all the relevant quantities will be obtained analytically, since the solution for the Lindblad master equation is available exactly. This model will serve as a motivating case to perform an in-depth analysis of the two main models presented in the manuscript: the quantum harmonic oscillator and a quantum Brownian particle. For the harmonic oscillator, since the Hilbert space characterizing the system has infinite dimensions, all the computations performed are potentially more complicated. For this reason, only some of the results are obtained analytically. Finally, the results presented for a quantum Brownian particle are obtained numerically.

V A simple case: Thermal qubit

Let us start with a preliminary analysis of the simplest system of interest: a two-level system, also called a qubit, weakly coupled to a thermal bath. Despite its simplicity, this is a paradigmatic example as the coupling of few level systems to thermal baths has been mastered in the last decades in different condensed matter platforms e.g., semiconductor quantum dots [66, 67] or superconducting qubits [63, 64, 65]. They are important pieces in the development of modern quantum thermodynamic engines [110, 13, 14, 15]. This simple case provides us with analytical understanding of the problem. It is important to remark that all the final conclusions drawn for the more involved examples of the harmonic oscillator and quantum Brownian motion will be in accordance with the ones obtained from this simple analysis.

Consider a two-level system weakly coupled to a thermal bath at inverse temperature $\beta$. Transitions between the ground ($n=0$) and the excited ($n=1$) states, split by an energy $\hbar\omega$, occur with rates $W_{10}=\gamma\bar{n}(\omega,T)$ and $W_{01}=\gamma[1+\bar{n}(\omega,T)]$ induced by the bath [96], with the coupling rate $\gamma$ and an average number of photons with frequency $\omega$ in a bath at temperature $T$, $n(\omega,T)$, given by the Bose-Einstein distribution

$$\bar{n}(\omega,T)=\left[\exp(\hbar\omega/k_{B}T)-1\right]^{-1}.$$ (21)

When thermalized, the state of the system can be written as a vector formed by the diagonal elements of the density matrix giving the occupation of the two states, $\rho=(\rho_{00}\ \rho_{11})^{T}$ in the Fock-Liouville representation. Then, the Lindblad equation, Eq. (2), is a simple rate equation

$$\displaystyle{\dot{\rho}(t)=\left(\begin{array}[]{cc}-\gamma\,\bar{n}(\omega,T)&\gamma\,[\bar{n}(\omega,T)+1]\\ \gamma\,\bar{n}(\omega,T)&\gamma\,[\bar{n}(\omega,T)+1]\end{array}\right)\rho(t).}$$ (24)

Note that, in the absence of coherence in the initial state, the Hamiltonian term of the Lindblad equation (2) does not contribute and the dynamics is purely dissipative.

We are interested in the relaxation from an initial thermal state at temperature $T_{0}=T+\Delta T=1/k_{\text{B}}\beta_{0}$. The time evolution of the density matrix can be obtained solving Eq. (24)

$$\rho(t) = \rho_{\beta}^{\rm th}+\frac{e^{-\Gamma t}(e^{\hbar\omega\beta}-e^{\hbar\omega\beta_{0}})}{(1+e^{\hbar\omega\beta})(1+e^{\hbar\omega\beta_{0}})}\left(\begin{array}[]{c}-1\\ 1\end{array}\right),$$ (27)

with the total rate $\Gamma=\gamma(1+2\bar{n}(\omega,T))=\gamma\coth(\hbar\omega\beta/2)$, that is is proportional to the thermal fluctuations of the bath. The fact that there is a single relaxation channel in this case emphasizes the role of the bath fluctuations (larger for higher temperatures). It is hence clear that the evolution will be faster when the system relaxes to a hotter steady state.

Figure 2: Thermalization kinematics for a qubit. (a) Dependence of the decaying mode corresponding to the non-zero eigenvalue. (b) Overlap of a state thermalized at a temperature $T+\Delta T$ with the stationary state at a temperature $T$, for different values of $\Delta T/T$. The black line in (b) corresponds to the asymptotic behaviour at large $\Delta T$.

Figure 3: Three-temperatures protocol for a qubit. (a) Bures distance as a function of $\Delta T$. The blue and red dots mark the equidistant temperatures considered for the cooling and heating protocols, respectively. These are used to compute (b) the Fidelity and (c) the velocity as function of time, and (d) the level of completion for a given time renormalized by $t_{\rm fin}=2/\gamma$. In all panels $\hbar\omega=1$, $T_{W}=\hbar\omega/2k_{\text{B}}$, and $\Delta T=T_{H}-T_{W}=\hbar\omega/k_{\text{B}}$, with $T_{C}\approx 0.70\hbar\omega$ chosen to be equidistant from $T_{W}$ as depicted in (a). For the time evolution $t_{0}=0$ in all cases. The dashed line in (a) is a quadratic expansion of the Bures distance around $\Delta T=0$, according to Eq. (56).

It is however convenient to look further into the details of the dynamics, as introduced in Sec. II, since it allows for an analytical treatment. We start by obtaining the eigenvalues of $\mathcal{L}$. In this simple case, the spectrum is reduced to only two values: $\lambda_{1}=0$ related to the trivial stationary state, and $\lambda_{2}=-\Gamma$, which takes into account the decay mode and depends on the bath parameters, contained in $\bar{n}(\omega,T)$ and in the coupling $\gamma$. The temperature dependence of the decaying mode corresponding to $\lambda_{2}$ is plotted in Fig. 2(a). Their corresponding (right) eigenvectors, see Eq. (5), are given by

$$\displaystyle{\Lambda_{1}^{r}=\frac{1}{e^{\beta\hbar\omega}+1}\left(\begin{array}[]{c}e^{\beta\hbar\omega}\\ 1\end{array}\right),}$$ (30)

corresponding to the stationary state of the system, characterized by $\lambda_{1}=0$, and

$$\displaystyle{\Lambda_{2}^{r}=\left(\begin{array}[]{c}1\\ -1\end{array}\right),}$$ (33)

which is related to the decaying mode, $\lambda_{2}=-\Gamma$. These results agree with the exact evolution obtained in Eq. (27).

To analyze the dynamics and the influence of the initial state, we need to compute the overlap between $\Lambda_{2}^{\ell}$ and the initial state, see Eq. (8). As we are starting from thermal equilibrium, the initial state coincides with the stationary state at temperature $T_{0}=T+\Delta T$, meaning that $\rho_{0}=\Lambda_{1}^{r}(T+\Delta T)$. Then, using

$$\Lambda_{2}^{\ell}=\frac{1}{1+e^{\hbar\omega\beta}}(1\ -e^{\hbar\omega\beta}),$$ (34)

the overlap is

$$\xi=\Tr(\Lambda_{2}^{\ell}\rho_{0}(T))=\frac{e^{\beta_{0}\hbar\omega}-e^{\beta\hbar\omega}}{\left(e^{\beta_{0}\hbar\omega}+1\right)\left(e^{\beta\hbar\omega}+1\right)}.$$ (35)

Note that the overlap has the same modulus under the exchange of temperatures $T$ and $T_{0}$. Hence, in a two-temperatures protocol, the overlap is the same in both ways (cooling and heating): $|\xi(T_{H},T_{C})|=|\xi(T_{C},T_{H})|$. As a consequence, any asymmetry in this protocol is to be attributed only to the monotonous increase of the coupling rate shown in Fig. 2(a): $\gamma$ is larger when relaxing to a hot bath, so $|\lambda_{2}^{H}|>|\lambda_{2}^{C}|$. As shown in Fig. 2(b), the overlap increases monotonically with $\Delta T$, i.e., far from equilibrium states are more strongly overlapped. Remarkably there is a suppression in $|\xi|$ for low temperature states that agree with Nernst’s unattainability principle [11]: the gap in the overlap at low $T$ and low $\Delta T$ prevents a cold system to be further cooled down, cf. the blue curve in Fig. 2(b). Note also that this Nernst gap is accompanied by a low-temperature plateau of the decaying rate, with $|\lambda_{2}|\to\gamma$, see Fig. 2(a). In the opposite case, the overlap of an infinite and a zero temperature states is maximal: $|\xi|\to 1/5$, with the bound $|\xi(T+\Delta T,T)|\leq|\xi_{\rm asym}|$ for the asymptotic value

$$|\xi_{\rm asym}|\to\frac{1}{2}\tanh\left(\frac{\hbar\omega}{2k_{B}T}\right)$$ (36)

when $\Delta T\gg T$, see Fig. 2(b).

To compute the fidelity as a measure of distance between two thermal states, let the system be described by a Gibbs state at an inverse temperature $\beta$ and frequency $\omega$ in the Fock-Liouville space, $\rho^{\text{th}}_{\beta}$ given by Eq. (30). Consider two thermal states at different inverse temperatures $\beta_{1}$ and $\beta_{2}$. In this case, the fidelity is simply

$$F(\rho_{\beta_{1}}^{\rm th},\rho_{\beta_{2}}^{\rm th})=\frac{1+e^{\hbar\omega(\beta_{1}+\beta_{2})/2}}{[(1+e^{\beta_{1}\hbar\omega})(1+e^{\beta_{2}\hbar\omega})]^{1/2}}.$$ (37)

With this expression, we calculate the Bures distance $D^{2}_{B}=2[1-F(\rho_{\beta_{1}}^{\rm th},\rho_{\beta_{2}}^{\rm th})]$, see Eq. (13), which we plot in Fig. 3(a). In Sec. VI.3 a similar analysis of the Bures distance and its dependence with the temperature difference is made for the harmonic oscillator and quantum Brownian particle cases.

As the density matrix time evolution is given by Eq. (27), we can find an analytical expression for the time evolution of the fidelity of a system initially at a state $\rho_{\beta_{0}}^{\rm th}$ with respect to the stationary state as it is put in contact with a bath at inverse temperature $\beta$, namely

$$F[\rho(t),\rho_{\beta}^{\rm th}]=\frac{\sqrt{e^{\hbar\omega\beta}[e^{\hbar\omega\beta_{0}}+A_{\beta\beta_{0}}(t)]}+\sqrt{1-A_{\beta\beta_{0}}(t)}}{(1+e^{\hbar\omega\beta})(1+e^{\hbar\omega\beta_{0}})},$$ (38)

where the time dependence is encapsulated in the term $A_{\beta\beta_{0}}(t)\equiv{(1-e^{-\Gamma t}})(e^{\beta\hbar\omega}-e^{\beta_{0}\hbar\omega})/(1+e^{\beta\hbar\omega})$. Note that at time $t=0$, $A(0)=0$, thereby recovering the fidelity given by Eq. (37). The fidelity as a function of time is plotted in Fig. 3(b) for the heating up and cooling down processes. This result confirming that the heating protocol is faster than the cooling one. We have verified that this effect is much stronger for low temperatures.

Even in this simple case the dynamics are less trivial than what one may suspect from the above discussion. This is revealed by using the quantum Fisher information to compute the thermal kinematic distance and speed, see Eqs. (16) and (17). In this case, being the density matrix diagonal, and $\partial_{t}\rho_{00}(t)=\partial_{t}\rho_{11}(t)$, the QFI reads ${\cal I}_{Q}[\rho(t)]=[\partial_{t}\rho_{00}(t)]^{2}/\rho_{00}(t)\rho_{11}(t)$, where we have also used $\rho_{00}(t)+\rho_{11}(t)=1$, leading to

$${\mathcal{I}}_{Q}[\rho(t)]=\frac{\Gamma^{2}}{[e^{\hbar\omega\beta}\kappa(t)-1][\kappa(t)+1]},$$ (39)

where the time-dependence is encapsulated in the term

$$\kappa(t)\equiv\frac{1+e^{\hbar\omega\beta_{0}}}{e^{\hbar\omega\beta}-e^{\hbar\omega\beta_{0}}}e^{\Gamma t},$$ (40)

and the instantaneous statistical velocity is

$$v(t)=\frac{\Gamma/2}{\sqrt{e^{\hbar\omega\beta}\kappa(t)-1}\sqrt{\kappa(t)+1}}.$$ (41)

This expression is displayed in Fig. 3(c). The velocity of the process is initially larger for the heating mechanism, confirming our hypothesis, but eventually it becomes slower than the cooling one. It can be concluded that as the system gets closer to thermalization the dynamics get slower. To analyze the full process, we also compute the statistical length

$$\begin{aligned} \ell(t_{0},t)&=\frac{1}{\Gamma}\left\{\arctan\left[\left|\frac{(e^{\hbar\omega\beta}{-}1)\kappa(t)-2}{2\sqrt{[e^{\hbar\omega\beta}\kappa(t){-}1][\kappa(t){+}1]}}\right|\right]\right.\\ &-\left.\arctan\left[\left|\frac{(e^{\hbar\omega\beta}-1)\kappa(t_{0})-2}{2\sqrt{[e^{\hbar\omega\beta}\kappa(t_{0}){-}1][\kappa(t_{0}){+}1]}}\right|\right]\right\}.\end{aligned}$$ (42)

with which we plot the ratio $\varphi(t)=\ell(0,t)/\ell(0,t_{\rm fin})$ in Fig. 3(d). Due to the initial higher velocity for the heating process, the degree of completion is larger at all times for the heating protocol.

Figure 4: Simulation of the protocols for the harmonic oscillator. (a) Bures distance as a function to the temperature difference with respect to the equilibrium state at $T_{W}$. The blue and red dots mark the initial temperatures $T_{C}$ and $T_{H}$ for the cooling and heating protocols, respectively. (b) Fidelity for the state of the system with respect to the thermal state at $T_{W}$. The process corresponds to the three-temperature protocol, evolving from $T_{C}$ to $T_{W}$ (orange line) and $T_{H}$ to $T_{W}$ (blue line). (c) Instantaneous statistical velocity for the harmonic oscillator computed for the heating and cooling protocol in the two-temperature scheme. (d) Degree of completion for the harmonic oscillator in the two-temperature scheme. In all plots, the temperature ranges are such that $n(\omega,T_{C})=1$ and $n(\omega,T_{H})=10$, the coupling parameter is $\gamma=0.1$, $\hbar\omega=1$, and $t_{\textrm{fin}}=50$.

VI Increasing complexity

VI.1 Quantum harmonic oscillator

After having introduced the main concepts presented in the manuscript with a clear and analytically solvable case, we will perform a similar analysis for more complicated and richer systems. The quantum harmonic oscillator allows us to derive some analytical expressions for the behavior of the system but, on the other hand, numerical methods are required to compute quantities of interest such as the quantum speed through the quantum Fisher information.

A harmonic oscillator is described by the following Hamiltonian

$$H = \hbar\omega a^{\dagger}a,$$ (43)

where $\hbar$ is the Planck’s constant, $\omega$ is the oscillator frequency and $a,a^{\dagger}$ are the annihilation and creation bosonic operators, respectively. The interaction with a thermal bath is described by the jump operators

$$\begin{aligned} L_{+}&=\sqrt{\gamma\,\bar{n}(\omega,T)}\;a^{\dagger},\\ L_{-}&=\sqrt{\gamma\,[\bar{n}(\omega,T)+1]}\;a,\end{aligned}$$ (44)

where $\gamma$ is the coupling strength with the bath, and $n(\omega,T)$ is the average number of excitations in the bath at a given temperature $T$ as in the two-level system case, see Eq. (21). The state of a thermal harmonic oscillator is simply determined by its average number of excitations, $n(\omega,T)$, as

$$\begin{aligned} \rho^{\text{th}}_{\beta}&=\sum_{n=0}^{\infty}\dfrac{\left[\bar{n}((\omega,T)\right]^{n}}{[1+\bar{n}((\omega,T)]^{n+1}}\ket{n}\bra{n}\\ &=\sum_{n=0}^{\infty}e^{-n\omega\beta}\left(1-e^{-\omega\beta}\right)\ket{n}\bra{n},\end{aligned}$$ (45)

being $\ket{n}$ the pure state of a system with $n$ photons. Note that as $\bar{n}(\omega,T)$ depends on the temperature and the frequency in terms, in the thermal state this dependence is implicitly included. Initially, we shall consider the system to be in such thermal state.

If a system is in a Gaussian state, including a thermal state, and the interaction with the bath is also Gaussian, its state would be entirely characterized by the evolution of its occupation numbers $\left<a\right>$ and $\left\langle a^{\dagger}a\right\rangle$. This fact reduces the problem to the computation of the evolution of the expected values instead of the whole density matrix, leading to a single ordinary differential equation. The dynamics of $\langle a\rangle$ and $\left\langle a^{\dagger}a\right\rangle$ are described by the following expressions [111]

$$\frac{d\langle a\rangle}{dt} = -i\left(\omega+\dfrac{\Gamma}{2}\right)\langle a\rangle,$$

$$\frac{d\left\langle a^{\dagger}a\right\rangle}{dt} = -\Gamma\left\langle a^{\dagger}a\right\rangle+\Gamma\bar{n}(\omega,T).$$ (46)

[being $\bar{n}(\omega,T)$ the average number of excitations of the bath in resonance with the oscillator. Focusing on the temporal evolution of the average number of the system excitations, the solution to this differential equation can be obtained, leading us to the variation in the average number of excitations

$$\left<a^{\dagger}a\right>_{t}= \left<a^{\dagger}a\right>_{0}e^{-\Gamma t}+\int_{0}^{t}\Gamma n(\omega,T)e^{-\Gamma(t-s)}ds,$$ (47)

Note that the sub-index $t$ represents the time dependence and $0$ the initial value for the average number of excitations. We consider that the system and the bath are in contact at $t=0$, without loss of generality. The protocol is modeled by a quench, i.e. a step function, so the $n(\omega,T)$ term in the integral is constant. Therefore, the evolution reduces to

$$\left<a^{\dagger}a\right>_{t}=\left<a^{\dagger}a\right>_{0}e^{-\Gamma t}+n(\omega,T)\left(1-e^{-\Gamma t}\right).$$ (48)

This expression for $\left<a^{\dagger}a\right>_{t}$ directly provides us with the temporal evolution of the average number of excitations of the harmonic oscillator. To compare the state of the system with the final thermal state, whose average number of excitations is the one of the bath, we analytically calculate the fidelity for the system at a time $t$, starting at time $t=0$ with an occupation number $\left<a^{\dagger}a\right>_{0}$ compared to a final thermal state at temperature $T$, defined by $\bar{n}(\omega,T)$. Then, the Fidelity reads

$$F(\rho^{\text{th}}_{S}(t),\rho_{\text{F}})=\Tr\sum_{n=0}^{\infty}\left[\dfrac{\left(\left<a^{\dagger}a\right>_{t}\bar{n}(\omega,T)\right)^{n}}{(1+\left<a^{\dagger}a\right>_{t})^{1+n}(1+\bar{n}(\omega,T))^{1+n}}\right]^{1/2}\ket{n}\bra{n} =\dfrac{1}{[(1+\left<a^{\dagger}a\right>_{t})(1+\bar{n}(\omega,T))]^{1/2}}\dfrac{1}{1-r}$$ (49)

where $r=\left\{\left<a^{\dagger}a\right>_{t}\bar{n}(\omega,T)/[(1+\left<a^{\dagger}a\right>_{t})(1+\bar{n}(\omega,T))]\right\}^{1/2}$. In Fig. 4 both the Bures distance as a function of the temperature, panel 4(a), and the Fidelity as a function of time, panel 4(b), are displayed. The results confirm behavior obtained for the two-level system. One main difference is that to obtain the same distance, in the Bures sense, we need a higher temperature difference in the harmonic oscillator case that for the two-level system. This is due to the infinite size of the Hilbert space of the harmonic oscillator, in comparison to a two-dimensional Hilbert space.

To analyze the thermal kinematics of the system we have computed numerically the quantum Fisher Information of the two-temperatures protocol. The results are displayed in Figure 4, panels 4(c) and 4(d). It is clear that, even if the harmonic oscillator is a different and more complicated system its thermal behavior is similar to the one for the two-level system. This supports our hypothesis that the thermal asymmetry is a general trend of open quantum systems. In the next section we check this behavior with an even more complex system as the Quantum Brownian particle.

Figure 5: Simulation of the protocols for the Quantum Brownian particle. (a) Bures distance dependence with the temperature difference with respect to the equilibrium state, as in the harmonic oscillator case. (b) Fidelity for the state with respect to the thermal state at $T_{W}$. The process corresponds to the three-temperature protocol, evolving from $T_{C}$ to $T_{W}$ (orange line) and $T_{H}$ to $T_{W}$ (blue line). (c) Instantaneous statistical velocity for the quantum Brownian particle computed for the heating and cooling protocol in the two-temperature scheme. (d) Degree of completion for the quantum Brownian particle in the two-temperature scheme. In the case of the quantum Brownian particle, the values of the parameters are $t_{\textrm{fin}}=5\cdot 10^{3},m=1,\Omega=10^{-3},\Lambda=1,\zeta=0.1$.

Figure 6: Near-temperature simulations for the three-temperature protocol for the harmonic oscillator. The Fidelity is calculated between the time evolution of the state of the system, $\rho(t)$, and the thermal state at $n(\omega,T_{W})$. In all cases $\omega=1$, the average excitation number for the cold bath is $n(\omega,T_{C})=1$ in all cases, and the hot temperature is for panels (a), (d) $n(\omega,T_{H})=1.1$; for panels (b), (e) $n(\omega,T_{H})=2$; and for panels (c), (f) $n(\omega,T_{H})=5$. In the linear regime, for close temperatures $T_{C}$ and $T_{H}$, both curves collapse and the processes remain symmetric. However, as the temperature difference increases, the lines of heating up and cooling down evolve differently creating the asymmetry between the protocols.

VI.2 Quantum Brownian particle in a trap

The third model introduced to analyze the protocols is a quantum Brownian particle trapped in a harmonic trap, following the results and experiments already performed in the classical case [54]. This is the most sophisticated case that is treated in the article, where all the relevant quantities need to be computed numerically.

A quantum Brownian particle interacting with a bosonic bath is described by the following Hamiltonian [112, 113, 89, 90]

$$H = H_{S}+H_{B}+H_{I}$$ (50)

$$= \frac{p_{S}^{2}}{2m_{S}}+\sum_{i=1}^{n}\frac{\kappa^{2}_{i}}{2m_{B_{i}}\omega_{B_{i}}}x_{S}^{2}+\phi(x_{S})$$

$$+ \sum_{i=1}^{n}\hbar\omega_{B_{i}}a^{\dagger}_{B,i}a_{B,i}-\sum_{i=1}^{n}\kappa_{i}x_{B,i}x_{S},$$

where the indexes $S,B$ hold for the system and bath operators respectively. Here $m_{S}$ is the mass of the Brownian particle, $x_{S}$ its position, $p_{S}$ its momentum and $\phi(x)$ is a trapping potential. Similarly, $m_{B_{i}}$, $\omega_{B_{i}}$ and $x_{B_{i}}$ are the mass, frequency, and position of the $i$-th bath particle, for all $i=1,\ldots,n$. The factors $\kappa_{i}$ represent the coupling between the system and the $i$-th bath mode.

The trapping potential is customarily taken as a harmonic term so that

$$\phi(x)=\dfrac{m}{2}\tilde{\omega}^{2}x^{2},$$ (51)

for a given trap frequency $\tilde{\omega}$.

A treatment for the problem can be performed by a Lindblad-equation-like transformation of the equations of motion of the system. The global evolution of the system and bath may be described by a unitary operator, and the state of the system at a given time $t$ is described by

$$\dot{\rho}_{S}(t)=\Tr_{B}\left\{U(t)\left(\rho_{S}(0){\otimes}\rho_{B}\right)U^{\dagger}(t)\right\}\equiv\mathcal{L}[\rho_{S}(t)],$$ (52)

being $\mathcal{L}$ the Liouvillian superoperator of the coherent dynamics. Given the fact that the interaction between the system and the bath is linear and assuming it to be also weak, we can consider a single Lindblad operator $\mathcal{A}$ such that [113]

$$\mathcal{A}(T) = \alpha x+\beta p,$$ (53)

for some parameters $\alpha,\beta\in\mathbb{C}$. In order to match the coefficients represented in Eq. (53) with a general Born-Markov treatment of the problem within the Caldeira-Leggett limit, these $\alpha$ and $\beta$ must be given by

$$\alpha=\dfrac{\left(2m\zeta k_{B}T\right)^{1/2}}{\hbar},$$ (54)

and

$$\operatorname{Re}(\beta)=-\dfrac{\zeta k_{B}T}{\hbar^{2}\alpha\Lambda},\quad\operatorname{Im}(\beta)=\dfrac{\zeta}{2\hbar\alpha}.$$ (55)

In these relations, $\Lambda$ is the so-called Lorentz-Drude cutoff appearing in baths with Ohmic spectral density; $\zeta$ is a damping constant, whose inverse is related to the relaxation scales; $T$ is the temperature of the bath and $m$ the mass of the oscillators. The Caldeira-Leggett limit is satisfied for large temperature and cut-off limits. Under this regime, one recovers the Caldeira-Leggett equation for general diffusion processes in a quantum framework for a quantum Brownian particle [89, 90]. Note that the temperature of the baths is related to the average number of excitations via the Bose-Einstein relation as in the previous models.

For this model both the fidelity and quantum Fisher information may be calculated numerically, as an analytical description is intractable. In Fig. 6 we observe a similar analysis to the ones performed for both the two-level system and the harmonic oscillator. In this case, the temperature range that we need to consider is even larger, due to the complexity of the bath. All the results are similar, confirming the general character of our results. One interesting feature is that during the heating up process in the three-temperatures protocol [c.f. Fig. 6(b)] the fidelity reaches a value close to one in a finite time, and then bounces down. This interesting behavior suggest that the system suffers from hysteresis, an interesting feature specially due to the Markovian character of the dynamics.

VI.3 Analysis of the results

In the three-temperature protocol, the fidelity of the state at a given time $t$ has been compared to the thermal state at the intermediate temperature, $T_{W}$, so that it increases to one, when thermalization takes place. As it is indicated in IV.1, the thermal state at warm temperature, $\rho^{\text{th}}_{\beta_{W}}$, is chosen to be equidistant to the cold, $\rho^{\text{th}}_{\beta_{C}}$, and hot states, $\rho^{\text{th}}_{\beta_{H}}$, so that the Bures distance is the same in both cases. The distance is depicted in Figs. 3(a), 4(a), and 5(a) as a function of the temperature difference $\Delta T$ with respect to $T_{W}$ , meaning that $\Delta T_{C}=T_{C}-T_{W}$ for the initial cold state (heating-up protocol), and $\Delta T_{H}=T_{H}-T_{W}$ for the initial hot state (cooling-down protocol). The temperature of both initial points is represented by a hollow red circle for the case of heating and by a hollow blue circle for the cooling. Both points evolve to the equilibrium thermal state, clearly represented by a minimum. It is worth noticing that the asymmetry in the different protocols is clearly appreciated here. Even if they are starting from the same initial point, the length of the path followed to the equilibrium state is clearly larger in the cooling protocol.

The asymmetry in the three-temperature protocols is analyzed by the use of the fidelity between the initial states and the target one, as a function of time. This is displayed in Figs. 3(b), 4(b), and 5(b). Remarkably, the necessary times to thermalize differ by several orders of magnitude due to the increasing complexity of each configuration. However, the three of them show a similar behavior. This fact is also relevant for the velocity analysis.

Regarding the two-temperature protocol, we shall make use of the instantaneous velocity quantity, Eq. (17), and degree of completion, Eq. (19). Figures 3(c), 4(c) and 5(c) represent this quantity for the qubit, the harmonic oscillator, and the Brownian particle respectively. As expected, the velocity in both cases for the heating process $(T_{C}\rightarrow T_{H})$ is larger than the cooling counterpart until $(T_{H}\rightarrow T_{C})$ it reaches a null value, i.e. the system thermalizes at the final temperature. As we anticipated in the previous paragraph, the thermalization times differ by several orders of magnitude. This is also represented by the scale in the velocity axis in both of the aforementioned plots. Finally, the computation of the degree of completion is directly derived considering the values of the instantaneous velocity, Eqs. (18) and (19). Figures 3(d), 4(d) and 5(d) show the temporal evolution of the degree of completion and, as it is expected, the functions are similar regardless of the sort of system, and the heating process takes less amount of time to be completed than the cooling. The main implication in the function is related to the quicker saturation to the unity of the function in the former case.

Although all methods based on thermal kinematics theory are useful for analyzing the dynamics towards the equilibrium state, they do not provide theoretical insight into its origin. The Bures distance with respect to the equilibrium temperature in the three-temperature case (see Figs. 3(a), 4(a), and 5(a)) offers intuition about the protocols. However, it is solely a measure of the distance and does not offer any fundamental reasoning; it is merely a consequence of a more fundamental phenomenon.

Figure 7: Eigenvalues of the Liouvillian operator $\mathcal{L}$ for the harmonic oscillator (a), and a quantum Brownian particle (c) for $\bar{n}(\omega,T_{H})=10$ (red points) and $\bar{n}(\omega,T_{C})=1$ (blue points). Note the difference of scales in the $x$-axis. Panels (b) and (d) represent the first eigenvalues of the respective panels (a) and (c) with a size proportional to the overlap with the thermal state at the opposite temperature for the harmonic oscillator, see Eq. (8). In both cases, the truncated dimension of the Hilbert space is $N=150$.

VI.4 Linear response regime

We shall conclude the section providing a brief comment on the near-equilibrium, i.e. linear regime, for thermal evolution close to the equilibrium temperature in the three-level protocol.

The linear response theory, developed mainly by Kubo [5, 6], is the cornerstone to analyze the near-equilibrium behavior in non-equilibrium classical thermodynamics. The fluctuation-dissipation theorem states that the fluctuation properties of a system in TE determine its linear response to an external perturbation [6]. In the quantum counterpart, this theorem has been derived for closed quantum systems and recently for open quantum systems [114, 115]. This extension allows us to apply the existing results from isolated equilibrium systems to open systems, with Lindbladian dynamics [114, 115]. Within this regime, one expects to recover the same as in classical thermodynamics results, where the asymmetry between heating and protocols is absent. That means, for a small temperature differences in both protocols, we expect to observe how the asymmetry diminish, since we are approaching the equilibrium state.

For the qubit case this phenomenon is clearly appreciated in the analytical derivation of the Fidelity comparing two states, Eq. (37). For small $\Delta T$, the Fidelity is quadratic

$$F(\rho_{\beta_{0}}^{\rm th},\rho_{\beta}^{\rm th})=1-\frac{e^{\beta\hbar\omega}}{8\left(1+e^{\beta\hbar\omega}\right)^{2}}\left(\frac{\hbar\omega\Delta T}{T^{2}}\right)^{2}+{\cal O}\left(\frac{\Delta T}{T}\right)^{3},$$ (56)

i.e. no asymmetry is expected for states close to the equilibrium.

Regarding the simulations for the harmonic oscillator, the results for close temperatures is depicted in Fig. 6(a) and 6(d). As the temperature difference increases, the asymmetry starts to appear, making this discrepancy in both protocols more acute the larger is this gap. Figures 6(b)-(c) and 6(e)-(f) show this behavior. Figures 6(a)-(c) represent the fidelity with respect to the thermal state at warm temperature $T_{W}$, i.e. in the three-temperature protocol. Similarly, Figs. 6(d)-(f) showcase the asymmetry in the two-temperature scenario, displayed in the velocity needed to reach the opposite state. It is clear that the asymmetry arises as one deviates from equilibrium when the temperature difference increases.

VII Spectral analysis

Performing a spectral analysis of the Liouvillian, we can gain intuition about the relaxation time of our models for the heating and cooling protocols, as the dynamics of the system is given by Eq. (7). In Fig. 7 the eigenvalues of both the harmonic oscillator and the Brownian particle are displayed. Due to the infinite size of the Hilbert space of the systems, we have used a truncated Fock basis of dimension $N$, large enough to display the general behavior. As discussed in Sec. II, the spectrum of the Lindbladian is composed of eigenvalues whose real part is negative, apart from the null eigenvalue which determines the stationary state. This decomposition does not depend on the initial state we consider but on the parameters defined in the Lindbladian, i.e. the Hamiltonian and jump operators, as well as the constants and variables defined therein. The fact that an open quantum system, when initialized at a given state, say $\rho_{0}$, evolves differently than from another initial state, say $\tilde{\rho}_{0}$, for the same Lindbladian depends on the overlap between the considered initial state and the left eigenvectors, see Eq. (8).

The spectra for the heating and cooling cases and the different systems are displayed in Fig. 7. For the Lindbladians with higher temperature there is a spreading of the eigenvalues towards the negative real axis, see Fig. 7(a) for the harmonic oscillator and Fig. 7(c) for the quantum Brownian particle. This means that, in the heating-up processes, there are many more fast-decaying modes than in the cooling-down counterparts, indicating that heating will be faster. This behavior is in agreement with Fig. 2(a) for the thermal qubit. In Figs. 7(b) and 7(d), the slowest decay modes, which act as bottlenecks to the dynamics, are plotted with the symbol size being proportional to the overlap with the initial state of the cooling down and heating protocols. The initial state in each case is chosen to be equal to a thermal state at the same temperature as the opposite process. This means that if the cooling-down/heating-up process is causing the system to evolve to $\bar{n}(\omega,T_{C})/\bar{n}(\omega,T_{H})$, the initial state will be $\rho^{\text{th}}_{\beta_{H}}/\rho^{\text{th}}_{\beta_{C}}$.

It is clear that in the cooling protocol, there is a higher overlap with slower decay modes in both cases. Moreover, the number of slower modes in the cooling case (represented by the blue stars) is larger than the number of modes in the heating case (red dots), apart from being closer between them and to the null eigenvalue. This spectral analysis provides us with a justification for the asymmetry in all the processes, referring all of them to mere observations of the decaying modes appearing in the spectra of the Liouvillians. For instance, the asymmetry appears naturally when computing the overlap of the state of the considered system with the generator of the dynamics, i.e. the Liouvillian: the ovelap is larger for the heating up protocol (blue stars) than for the cooling down (red dots). This explanation allows us to justify and clarify all the results obtained throughout the article. We recall however that this asymmetry in the overlap was not present in the thermal qubit case.

VIII Conclusions

In this paper, we have studied and unveiled an intriguing effect of non-equilibrium open quantum systems: the asymmetry of the time-evolution of heating up and cooling down trajectories. By introducing quantum information measures such as the fidelity, the Bures distance, and the quantum Fisher information we analyzed this phenomenon in two different protocols. The first protocol involves an intermediate temperature, equidistant between a hotter and a colder one, while the second protocol works between two absolute temperatures. The measures developed in this work are general and applicable to various other dissipative processes.

We extended the thermal kinematics to open quantum systems and applied these protocols to three different configurations of increasing complexity: a thermal qubit, a harmonic oscillator coupled to a bosonic heat bath, as well as a canonical model for the quantum Brownian motion. The qubit system provides an analytical description that can be solved exactly for all the studied magnitudes, suggesting a connection to the third law of thermodynamics; the other systems are analyzed numerically. Our results unequivocally indicate that heating up and cooling down are intrinsically different processes, with heating up always being the fastest. In the limit of small temperature differences we recover a symmetric behavior in accordance with equilibrium thermodynamics in the quantum regime.

By studying the Liouvillian spectrum of the system, we observe that the eigenvalues spread towards the negative real line as temperature increases. This indicates that for thermal baths at higher temperatures there are more fast decaying modes, making the evolution faster. Additionally, the overlap between the initial state and the fast-decaying modes confirms that the heating up is always a faster mechanism than the heating up.

Despite their simplicity, the proposed configurations can be readily be tested experimentally is various platforms e.g., semiconductor qubits [67] or superconducting cavity quantum thermodynamic circuits [65]. As systems with higher complexity require longer times to thermalize, harmonic oscillators or quantum Brownian motors are ideal canditates to detect thermalization asymmetries.

Appendix A Classical and Quantum Fisher Information

In classical parameter estimation a canonical measure is the classical Fisher Information, $I(\theta)$, of a probability density $p(x,\theta)$, defined as

$$I_{\text{cl}}(\theta):=\int_{-\infty}^{\infty}\left(\dfrac{d\log p(x,\theta)}{d\theta}\right)^{2}p(x,\theta)dx.$$ (57)

The geometric interpretation of the Fisher information arises from defining a statistical line element, $ds$, such that $ds^{2}:=I(\theta)d\theta^{2}$. Therefore, the line element $ds$ can be regarded as a dimensionless distance between probability densities $p(x,\theta)$ and $p(x,\theta+d\theta)$.

Consider a quantum state, $\rho$, parametrized by an $n$-dimension vector $\vec{\theta}=\left(\theta_{1},\ldots,\theta_{n}\right)$, and denoted by $\rho(\vec{\theta})$. For an infinitesimal change in the parameters, one can relate the Bures distance to the Quantum Fisher Information Matrix (QFIM), $\mathcal{I}$, so that

$$\left[D_{B}\left(\rho(\vec{\theta}),\rho(\vec{\theta}+d\vec{\theta})\right)\right]^{2}=\dfrac{1}{4}\sum_{i,j}\mathcal{I}_{ij}dx_{i}dx_{j}+\mathcal{O}(dx^{4}).$$ (58)

The complete derivation can be found in [109]. The elements of the QFIM are given by

$$\mathcal{I}_{ij}=\Tr\left[L_{\theta_{i}}\rho({\vec{\theta}})L_{\theta_{j}}\right],$$ (59)

where $\left\{L_{\theta_{i}}\right\}_{i=1}^{n}$ are the symmetric logarithmic derivative (SLD) operators for the $k-$th parameter, implicitly defined as

$$\dfrac{\partial\rho(\vec{\theta})}{\partial{\theta_{k}}}:=\dfrac{L_{\theta_{k}}\rho(\vec{\theta})+\rho(\vec{\theta})L_{\theta_{k}}}{2}.$$ (60)

We are only interested in single-parameter estimation, in this case the QFIM reads

$$\mathcal{I}_{\theta}=\Tr\left[L_{\theta}^{2}\rho(\theta)\right].$$ (61)

We are intended to obtain an operational expression for the SLD. In the eigenbasis of $\rho(\theta)$, by means of the spectral theorem, the density matrix can be decomposed in terms of its eigenvalues and eigenvectors

$$\rho(\theta)=\sum_{i=1}^{n}\lambda_{i}(\theta)\ket{\lambda_{i}(\theta)}\bra{\lambda_{i}(\theta)}.$$ (62)

Hence, in the eigenbasis of the state $\rho(\theta)$, the SLD operator is simply given by

$$L_{\theta}=2\sum_{i,j}\dfrac{\left\langle\lambda_{i}(\theta)\middle|\dfrac{d\rho(\theta)}{d\theta}\middle|\lambda_{j}(\theta)\right\rangle}{\lambda_{i}(\theta)+\lambda_{j}(\theta)}\left|\lambda_{i}(\theta)\right\rangle\left\langle\lambda_{j}(\theta)\right|,$$ (63)

where $\left\{\ket{\lambda_{k}(\theta)}\right\}_{k=1}^{n}$ is the eigenbasis of $\rho(\theta)$ for $\lambda(\theta)_{i}+\lambda(\theta)_{j}\neq 0,\ \forall i,j=1,\ldots,n$. For our analysis, we use only the time as a parameter, giving

$$L_{t}=2\sum_{i,j}\dfrac{\left\langle\lambda_{i}(t)\middle|\dfrac{d\rho(t)}{dt}\middle|\lambda_{j}(t)\right\rangle}{\lambda_{i}(t)+\lambda_{j}(t)}\left|\lambda_{i}(t)\right\rangle\left\langle\lambda_{j}(t)\right|,$$ (64)

and the QFI

$$\mathcal{I}_{Q}=\Tr\left[L_{t}^{2}\rho(t)\right].$$ (65)