Toward relaxation asymmetry: Heating is faster than cooling

An asymmetry in thermal relaxation toward equilibrium has been uncovered for Langevin systems near stable minima [Phys. Rev. Lett. 125, 110602 (2020)]. It has been shown that, given the same degree of nonequilibrium of the initial distributions, relaxation from a lower temperature state (heating) is faster than that from a higher temperature state (cooling). In this study, we elucidate this relaxation asymmetry for discrete-state Markovian systems described by the master equation. We rigorously prove that heating is faster than cooling for arbitrary two-state systems, whereas for systems with more than two distinct energy levels, the relaxation asymmetry is no longer universal. Furthermore, for systems whose energy levels degenerate into two energy states, we find that there exist critical thresholds of the energy gap. Depending on the magnitude of the energy gap, heating can be faster or slower than cooling, irrespective of the transition rates between states. Our results clarify the relaxation asymmetry for discrete-state systems and reveal several hidden features inherent in thermal relaxation.


I. INTRODUCTION
Systems attached to thermal reservoirs will relax toward a stationary state. Such thermal relaxation processes are ubiquitous in nature and possess rich properties from both dynamic and thermodynamic perspectives. One of the counterintuitive behaviors is the Mpemba effect [1], where cooling a hot system is faster than cooling a cold system. Such nonmonotonic relaxation phenomena have been observed in various systems [2][3][4][5][6] and theoretically analyzed for microscopic dynamics [7][8][9]. In addition, it was found that cooling a system before heating it could lead to exponentially fast relaxation [10]. From the perspective of thermodynamics, thermal relaxation processes exhibit universal relations regarding irreversibility, which is quantified by irreversible entropy production [11]. Notably, it has been shown that irreversible entropy production during thermal relaxation is lower-bounded by information-theoretical [12][13][14] and geometrical [15] distances between the initial and final states in both classical and quantum regimes. These relations imply stronger inequalities than the conventional second law of thermodynamics and impose geometrical constraints on the possible relaxation path. Since thermal relaxation is important in condensed matter [16] and heat engines [17], deepening our understanding of thermal relaxation would benefit research in these areas.
Consider preparing a thermal state corresponding to a given temperature via thermal relaxation; i.e., attaching a system to a single reservoir and allowing it to relax toward equilibrium. In this setting, the relaxation time is a quantity of interest and can be approximately estimated * tanvu@rk.phys.keio.ac.jp; Current Address: Department of Physics, Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama 223-8522, Japan † hasegawa@biom.t.u-tokyo.ac.jp via the convergence rate of the system state toward the equilibrium state [18]. Given two identical systems initiated in thermal states, one at a lower and the other at a higher temperature than the given temperature, a natural question arises: which one relaxes faster? Recently, this question has been addressed by Lapolla and Godec [19] for continuous-state Langevin systems. By considering a pair of thermodynamically equidistant temperature quenches (which have the same nonequilibrium free energy difference), they unveiled an unforeseen asymmetry in thermal relaxation; i.e., relaxation from a lower temperature is faster than that from a higher temperature. Roughly speaking, it implies that heating up cold objects is faster than cooling down hot objects. This phenomenon has been proven for quenches of dynamics near stable minima; however, it is not universal for generic systems because counterexamples have been constructed using multi-well potentials [19]. Another recent study [20] has reported that heating can occur faster or slower than cooling even for anharmonic single-well potentials, and a crossover region emerges if the quenches are not too far from equilibrium. Nonetheless, the relaxation asymmetry may be universal for a specific class of multi-well potentials. It is well known that overdamped diffusion under a multi-well potential with sufficiently high barriers converges to Markov jump dynamics at long times. By projecting the Langevin dynamics onto the Markov jump process between basins, it has been shown that the general asymmetry is preserved in degenerate potentials with separated time scales [19].
Since the average energy of a thermal state increases with the temperature, the relaxation asymmetry allows us to say that, from the energetic perspective, uphill relaxation is faster than downhill relaxation, which is counterintuitive to an extent. Moreover, relaxation speed cannot be characterized solely by thermodynamic quantities such as dissipation or frenesy [21]. Therefore, it is highly nontrivial that free energy plays an essential role as a arXiv:2102.07429v4 [cond-mat.stat-mech] 16 Nov 2021 quantifier of nonequilibrium degree in equidistant temperature quenches.
In this study, we elucidate the relaxation asymmetry for discrete-state systems modeled by Markov jump processes, thus improving our understanding of thermal relaxation. First, we prove that heating is faster than cooling in an arbitrary two-state system, affirming that there is universality in the relaxation asymmetry. However, we find that it is not the case for generic systems with at least three distinct energy levels. By analytically constructing counterexamples, we demonstrate that heating can be faster or slower than cooling, depending on the transition rates. Nevertheless, restricted to a particular class of systems, some universal results on the relaxation asymmetry are obtained. We show that when the energy levels of the system are two-state degenerate, there exist two critical energy gap thresholds. Depending on whether the energy gap is larger or smaller than these thresholds, it can be concluded with certainty that heating is faster or slower than cooling. These theoretical results are numerically demonstrated using several discrete-state systems.

II. SETUP
We consider the thermal relaxation process of an open system with N states. The system is coupled to a thermal reservoir at the inverse temperature where k B is the Boltzmann constant. Owing to interaction with the thermal reservoir, stochastic transitions between states are induced. The dynamics of the system is governed by the master equation, where |p t := [p 1 (t), . . . , p N (t)] denotes the probability distribution of the system at time t; the matrix R = [R mn ] ∈ R N ×N is time-independent with R mn ≥ 0 denoting the transition rate from state n to state m ( = n), and m R mn = 0. Without loss of generality, we assume that E 1 ≥ · · · ≥ E N , where E n denotes the energy of state n. The transition rates satisfy the detailed balance condition, R mn e −β f En = R nm e −β f Em , which is a sufficient condition such that the system always relaxes to the thermal Gibbs state |π f after a sufficiently long time, irrespective of the initial state. Here, For m = n, the transition rate R mn can be expressed as where B mn = B nm are the barrier coefficients, and Γ is a positive constant. Now, let us formulate the problem. We consider relaxation that initiates from a thermal state |π i associated with the inverse temperature β i = (k B T i ) −1 . This can be regarded as a temperature quench T i → T f at time t = 0 − . Given a pair of cold and hot temperatures, T c and T h , satisfying T c < T f < T h , we investigate the relaxation speed depending on the quench direction, T i = T c ↑ T f (heating) and T i = T h ↓ T f (cooling). The degree of nonequilibrium (or free energy) of each initial state is the same, where the relative entropy between two distributions |p and |q is given by D(p q) := n p n ln(p n /q n ).
Note that the characters, c and h, are associated with the initial temperatures; therefore, c and h correspond to the heating and cooling processes, respectively. For convenience, we define D(p) := D(p π f ). We aim to answer the question of which quench direction fastens the relaxation. To this end, we first need to quantify the relaxation speed, which can be evaluated by the distance between the system state and the thermal state. Analogous to Refs. [7,19], the relative entropy is used to measure the distance between states. In thermal relaxation, the relative entropy is closely related to the free energy and irreversible entropy production [22] as where F(p t ) = n p n (t)[E n + k B T f ln p n (t)] is the free energy of the distribution |p t , F i := −k B T i ln Z βi , and Σ t is the irreversible entropy production. Note that F(p 0 + ) = F i since F(p 0 + ) and F i are evaluated with temperatures T f and T i , respectively. Let |p i t be the time-evolution distribution corresponding to the initial state |π i , then heating is said to be faster (slower) than cooling if in long times. Throughout this study, we compare the relaxation speeds of the two quenches in the long-time regime, not over the entire evolution time. Therefore, it is possible that the heating and cooling curves, D(p c t ) and D(p h t ), may intersect at some point. In what follows, we explain in detail how to determine the relaxation speed in the long-time regime.
Let 0 = λ 1 > λ 2 > · · · > λ N be the eigenvalues of the transition matrix and {|v n } n be the set of corresponding eigenvectors, R|v n = λ n |v n .
Notably, all eigenvalues are real numbers since matrix R satisfies the detailed balance condition [23]. The eigenvectors {|v n } form a basis for the space R N with |v 1 = |π f , and |v n is a traceless vector for all n ≥ 2 (i.e., 1|v n = 0 since 1|R = 0|). Here, |0 (|1 ) denotes the N -dimensional vector with all zero (one) elements. Therefore, the initial distribution |π i can be expressed as a linear combination of {|v n } as follows: where γ i n 's are real numbers. Consequently, the probability distribution at time t can be analytically written in the following form: In the long-time limit, the probability distribution |p i t can be approximated up to the second-order term as Thus, the relaxation speed can be quantified via the value of |γ i 2 | [7]. Accordingly, heating is faster (slower) than cooling if |γ c 2 | < (>)|γ h 2 | (see Appendix A for proof). A closed form of γ i 2 can be obtained analytically [8]. The transition rate matrix R can be transformed to a symmetric matrix R = [R mn ] ∈ R N ×N as follows: The elements of matrix R can be explicitly written in terms of the elements of matrix R as Notably, matrix R has the same eigenvalues as R, and its eigenvectors {|f n } are related to those of R as |f n = F 1/2 |v n . Moreover, these eigenvectors are mutually orthogonal, f m |f n = v n |F|v m δ mn . Multiplying f 2 |F 1/2 on both sides of Eq. (10), we can show that γ i 2 is proportional to the inner product between the initial distribution |π i and the vector |f 2 , given by where |f 2 := F 1/2 |f 2 . Note that f 2 |π f = 0 since |f 2 and |f 1 = F 1/2 |π f are orthogonal. Because the sign of γ i n can be absorbed by changing the eigenvectors |v n → −|v n , hereafter, we assume γ h 2 ≤ 0, which implies that f 2 |π h ≤ 0.

III. RESULTS
Given the above setup, we now present our main results on the relaxation asymmetry, including numerical illustrations and proofs.
A. Two-state systems First, we consider two-state systems, for which universality regarding the relaxation asymmetry can be achieved.
Result 1. For two-state systems, heating is faster than cooling.
Result 1 affirmatively validates that the relaxation asymmetry is universal in two-state systems. Even for two-state systems, it is highly nontrivial that heating is faster than cooling. For discrete-state Markovian dynamics, the speed of state transformation is constrained by time-antisymmetric dissipation and time-symmetric frenesy (or dynamical activity) [24,25]. Relaxation from a thermal state at a higher temperature has a higher dynamical activity; more precisely, the average number of jumps over all the stochastic trajectories in the hot two-state system is always greater than that in the cold two-state system. Consequently, one may intuitively expect that cooling, which has higher dynamical activity, is faster than heating. However, the result is counterintuitive, implying that the dynamical activity alone cannot account for this relaxation asymmetry.
It is worth discussing the compatibility of Result 1 with the counterexamples in Ref. [19], in which cooling can be faster than heating for some double-well potentials. In the counterexamples therein, there is a time-scale separation: a local equilibration arises prior to the terminal exponential relaxation, which is effectively a two-state Markov jump process. Although heating occurs faster than cooling in the early stages, the asymmetry is inverted in the later stages. It has been observed that the breaking of the asymmetry is intimately related to the configuration between the intra-well and inter-well entropies [19]. Therefore, the apparent discrepancy between Result 1 and the terminal relaxation trend in the counterexamples is due to the intra-well entropic contribution, which is neglected in the Markov jump dynamics considered here.
To illustrate the above result, we use a two-state system [see Fig. 1(a)] with the following transition matrix: where ε = E 1 − E 2 > 0 is the energy gap. We vary the value of ε while fixing the inverse temperatures as follows: β f = 1, β h = 0.1, and β c is uniquely determined via the condition in Eq. (4). The time variation of ratio R t := D(p c t )/D(p h t ) is plotted as a function of time t in Fig. 1(b). Note that R t < (>)1 implies that heating occurs faster (slower) than cooling at time t. As shown, ratio R t is always smaller than 1, and there is no intersection between the heating and cooling curves at any finite time. Therefore, it is numerically verified that heating is always faster than cooling. Although Result 1 only indicates the long-time behavior, numerical evidence suggests that heating is faster than cooling for the entire evolution time.
Proof of Result 1. It suffices to prove that |γ c 2 | < |γ h 2 | = −γ h 2 . For N = 2, an arbitrary probability distribution |π i can be expressed as a point (x i , 1 − x i ) in the twodimensional space. Since E 1 > E 2 , all thermal states, |π c , |π f , and |π h , lie on the segment with (0, 1) and (1/2, 1/2) as endpoints. These probability distributions are geometrically illustrated in Fig. 1. In the following, we employ a geometrical approach to prove Result 1.
From the conditions, f 2 |π f = 0 and f 2 |π h ≤ 0, we can conclude that f 2 |π c ≥ 0 or γ c 2 ≥ 0. Thus, it is sufficient to show that which is equivalent to proving that (x c + x h )/2 > x f . We can rewrite the equality D(π c π f ) = D(π h π f ) as follows: where S(p) := − n p n ln p n is the Shannon entropy of the distribution |p . Note that x c < x h < 1/2 and is a strictly convex function over x i ∈ [0, 1/2] (i.e., the second derivative of g(x i ) with respect to x i is positive, g (x i ) > 0). Applying the Hermite-Hadamard inequality for g(x i ), we obtain the following: or equivalently, Combining Eqs. (18) and (21) results in the following inequality: Since g(x i ) is a strictly decreasing function, we have (x c + x h )/2 > x f , which completes the proof.

B. Systems with at least three distinct energy levels
Next, we consider more general systems that have at least three distinct energy levels; i.e., there exist three indices 1 ≤ i < j < k ≤ N such that E i > E j > E k . For such systems, we obtain the following result.

Result 2.
For systems with more than two distinct energy levels, heating can be faster or slower than cooling, depending on the transition rates.
Result 2 shows that no universality of the relaxation asymmetry is achieved in the general case. With an appropriate choice of barrier coefficients, we can construct a discrete-state system with |γ c 2 | < |γ h 2 | or |γ c 2 | > |γ h 2 | as desired. This difference between the continuousand discrete-state systems can be explained as follows. For simplicity, we consider a single-particle system described by the overdamped Langevin equation. In this continuous-state system, the particle tends to transit to a close place at any instant time. By contrast, in discretestate systems, the particle, in principle, can jump to anywhere, provided the transition rate between these states is positive. This degree of freedom could lead to a complicated relaxation compared with the continuous-state system. As shown below, the construction of the transition matrix that determines the magnitude relation between the heating and cooling rates is somewhat artificial; therefore, there may not be a realistic system with such a transition matrix. We anticipate that the universality of the relaxation asymmetry will be achieved if appropriate constraints are placed on the transition rates.
We numerically illustrate Result 2 in a three-state system [see Fig. 3(a)]. The transition rates are determined  Fig. 3(c), heating is faster than cooling in the shortand long-time regimes. However, interestingly, two crossing points exist in the intermediate-time regime, where cooling temporarily occurs faster than heating. While it may be difficult to identify these two crossing points in Fig. 3(b), we can clearly see in Fig. 3(c) that they appear around times t = 10 −1 and t = 20. In contrast, in the E 3 = 0.5 case, Fig. 3(e) shows that the R t ratio is always greater than 1, implying that cooling is faster than heating for the entire evolution time. These numerical results indicate that the relaxation asymmetry is not universal in three-state systems, which is consistent with the theoretical finding of Result 2. Here, for simplicity, we have performed numerical calculations with different energy levels in the two cases. Nevertheless, as shown below, even under the condition where the energy levels are fixed, it is easy to construct a transition matrix that yields the desired relaxation trend.
Proof of Result 2. We prove the result by analytically constructing a transition rate matrix such that heating is slower than cooling. A transition rate matrix for the opposite case can also be analogously constructed. First, one can prove that |π c , |π h , and |π f are linearly independent (see Appendix B). Let {|e 1 , |e 2 } be an orthogonal basis of the space S spanned by |π f and |π h , and P := |e 1 e 1 | + |e 2 e 2 | be the projection matrix to the space S. Define |f 2 := |π c − P|π c , then |f 2 = |0 because |π c , |π h , and |π f are linearly independent. Trivially, Next, we construct a transition rate matrix R that results in |γ c 2 | > |γ h 2 |. Set B mn = E m + E n , one can explicitly calculate that the matrix R has a single zero eigenvalue associated with the eigenvector |f 1 = F 1/2 |π f , and the remaining eigenvalues are all −Z β f [8]. Let U = {|x ∈ R N | x|f 1 = 0} be a subspace orthogonal to |f 1 . Then, there exists an orthogonal basis {|f 2 , |f 3 , . . . , |f N } of U, where |f 2 := F −1/2 |f 2 , since f 2 |f 1 = 0 (i.e., |f 2 ∈ U). Obviously, |f n (n ≥ 2) is an eigenvector of R with the corresponding eigenvalue −Z β f . Following the idea in Ref. [8], we slightly modify R as follows: Here, Z β f > 2 > · · · > N ≥ 0 are small numbers that ensure the positivity of R mn (m = n). It is easy to check that R|f 1 = |0 and R|f n = (−Z β f + n )|f n for all n ≥ 2. Now, the matrix R has N different eigenvalues, and |f 2 is precisely the eigenvector corresponding to the second-largest eigenvalue λ 2 = −Z β f + 2 . The transition rate matrix can be recovered as R = F −1/2 RF 1/2 , and the detailed balance condition is satisfied due to the symmetry of R. With this construction, the relation between |γ c 2 | and |γ h 2 | can be clarified as which completes the proof.

C. Degenerate two-level systems
Last, we consider the remaining case, wherein the energy levels are degenerate to two energy states. In other words, there exists an index 1 ≤ n < N such that E 1 = · · · = E n > E n+1 = · · · = E N . Such degenerated two-level systems are seen in atoms [26] and have been used to enhance quantum-annealing performance [27] and dissipation-less heat current [28]. For convenience, we define the energy gap ∆E := E n − E n+1 > 0. Remarkably, we find that, depending on the magnitude of this energy gap, heating can be faster or slower than cooling, regardless of the transition rates. Details are summarized in the following.
Note that Result 1 can be derived from Result 3 by setting N = 2 and n = 1. Result 3 indicates that there are two critical thresholds of the energy gap ∆E. As the energy gap is above or below these thresholds, a universal conclusion on asymmetry in thermal relaxation can be drawn. It is also highly nontrivial that the energy gap affects the relaxation speeds of heating and cooling in this way. When ∆E is large, the jump from energy state E n+1 to E n is less likely to occur compared with the opposite jump. Thus, heating is expected to be slower than cooling. However, counterintuitively, heating is faster than cooling as ∆E is sufficiently large. In addition, provided n ≤ N/2, heating is faster than cooling, regardless of the value of ∆E. This implies that the number of excited states also plays a crucial role in determining the relaxation speed.
Again, we numerically illustrate Result 3 in a threestate system with n = 2 [see Fig. 4(a)]. The transition rates are analogously defined as in Eq. (23) with γ 13 = γ 31 = 10, γ 23 = γ 32 = 0.1, and γ 12 = γ 21 = 0. We consider two cases of the parameter setting: (i) Here, κ ∈ (0, 1) is a tuning parameter. Note that once the above parameters are provided, the other parameters are uniquely determined. According to Result 3, cases (i) and (ii) correspond, respectively, to faster heating and faster cooling. We vary the value of κ from 0.1 to 0.5 and plot the numerical results of cases (i) and (ii) in Figs. 4(b) and 4(c), respectively. Figure 4(b) shows that R t < 1 in the long-time regime, implying that heating is faster than cooling. Interestingly, there are crossing points between the heating and cooling curves in the intermediate-time regime. In both the system considered here and the twostate system considered in the previous section, the energy levels are two-state degenerate. However, the former shows a more complicated relaxation trend than the latter. For case (ii), conversely, R t is always greater than 1, implying that cooling always occurs faster than heating. Consequently, all these numerical results are consistent with Result 3.
Proof of Result 3. We employ the same strategy used in proving Result 1. We prove the former case first, i.e., β h ∆E ≥ ln[n/(N − n)] leads to a faster heating. It can be observed that all points lying on the segment with |π c and |π h as endpoints are thermal states. It is also evident that |π c , |π f , and |π h are linearly dependent, i.e., there exists a real number a ∈ (0, 1) such that |π f = a|π c + (1 − a)|π h . Since f 2 |π f = 0 and f 2 |π h ≤ 0, f 2 |π c ≥ 0 or γ c 2 ≥ 0 follows from the continuity of the inner product f 2 |p for |p ∈ . Thus, it suffices to show that γ c 2 + γ h 2 < 0 or f 2 |(π c + π h )/2 < 0, which is equivalent to proving that (x c + x h )/2 > x f , where x i = π i n completely characterizes the thermal state |π i . The condition D(π c π f ) = D(π h π f ) can be rewritten as follows: Note that x c < x h ≤ 1/(2n) since β h ∆E ≥ ln[n/(N −n)], and g(x i ) := dS(π i )/dx i = n ln[(1 − nx i )/(N − n)x i ] is a strictly convex function over x i ∈ [0, 1/(2n)]. Applying the Hermite-Hadamard inequality for g(x i ) and following the same steps as in Eqs. (20) and (21), we obtain (x c + x h )/2 > x f , which proves the former case. When β c ∆E ≤ ln[n/(N − n)], one can derive that 1/(2n) ≤ x c < x h , and g(x i ) is a strictly concave function over x i ∈ [1/(2n), 1] (i.e., the second derivative of g(x i ) with respect to x i is negative, g (x i ) < 0). Applying the Hermite-Hadamard inequality for g(x i ), we obtain the following:

IV. SUMMARY AND DISCUSSION
In this study, we elucidated the relaxation asymmetry for discrete-state systems described by Markov jump processes. We proved that the relaxation asymmetry is universal in two-state systems, but not in generic systems with more than two distinct energy levels. For systems with two degenerate energy levels, we obtained some universal results indicating that the asymmetry in thermal relaxation depends on the energy gap and the number of excited states.
Notably, the relaxation asymmetry has recently been numerically studied for few-level open quantum systems described by the Lindblad master equations [29]. When the initial density matrix contains no coherence, the quantum systems can be described by classical Markov jump processes with the population distributions. It has been shown that heating is always faster than cooling for two-level systems, whereas heating can be faster or slower than cooling for three-and four-level systems. These numerical demonstrations are consistent with Results 1 and 2, thus affirmatively supporting our theoretical findings.
Although the relaxation asymmetry is not universal in generic discrete-state systems, universality may be achieved by imposing some constraints on the transition rates. Such a question requires further investigation and will be addressed in future work. It would also be interesting to study the relaxation asymmetry for non-Markovian systems [30]. For instance, the relaxation of some non-Markovian processes, such as a tagged particle in a single file and the end-to-end distance in a Rouse polymer, has been studied in Ref. [19]. Another possible direction involves formulating and investigating the relaxation asymmetry in open quantum systems [31,32] where coherence is present in the initial state. n=2 γ i n e λnt |v n . In the long-time limit, the term N n=2 γ i n e λnt |v n vanishes. Since D(p p + dp) = N k=1 dp 2 where ∆ = a mn e (λm+λn)t + O(e 3λ2t ), (A3) where a mn are constants. The first term on the righthand side is positive since |γ c 2 | < |γ h 2 |. The remaining terms may be negative; however, they are negligible compared with the first term in the long-time limit. Therefore, D(p c t ) < D(p h t ) as t → ∞.

Appendix B: Proof of linear independence
To prove the linear independence of |π c , |π f , and |π h , it is sufficient to show that the determinant of the following matrix is negative.
Here, E 1 > E 2 > E 3 and β c > β f > β h . Without loss of generality, we assume that E 3 = 0 and β h = 0. In this case, the determinant can be calculated as follows: (B2) Therefore, |X| < 0 is equivalent to Set h(x) := (1−e −βcx )/(1−e −β f x ), Eq. (B3) is equivalent to h(E 1 ) < h(E 2 ). We need only prove that h(x) is a strictly decreasing function over x > 0. Taking the derivative of h(x) with respect to x, we have we have dh(x)/dx < 0, which completes the proof.