The Mpemba index and anomalous relaxation

The Mpemba effect is a counter-intuitive relaxation phenomenon, where a system prepared at a hot temperature cools down faster than an identical system initiated at a cold temperature when both are quenched to an even colder bath. Such non-monotonic relaxations were observed in various systems, including water, magnetic alloys, polymers, and driven granular gases. We analyze the Mpemba effect in Markovian dynamics and discover that a stronger version of the effect often exists for a carefully chosen set of initial temperatures. In this \emph{strong Mpemba effect}, the relaxation time jumps to a smaller value leading to exponentially faster equilibration dynamics. The number of such special initial temperatures defines the \emph{Mpemba index}, whose parity is a topological property of the system. To demonstrate these concepts, we first analyze the different types of Mpemba relaxations in the mean field anti-ferromagnet Ising model, which demonstrates a surprisingly rich Mpemba phase diagram. Moreover, we show that the strong effect survives the thermodynamic limit and that it is tightly connected with thermal overshoot -- in the relaxation process, the temperature of the relaxing system can decay non-monotonically as a function of time. Using the parity of the Mpemba index, we then study the occurrence of the strong Mpemba effect in a large class of thermal quench processes and show that it happens with non-zero probability even in the thermodynamic limit. This is done by introducing the \emph{isotropic} model for which we obtain analytical lower bound estimates for the probability of the strong Mpemba effects. Consequently, we expect that such exponentially faster relaxations can be observed experimentally in a wide variety of systems.


I. INTRODUCTION
The physics of thermal relaxation is rich with fascinating and often surprising behaviors. A particularly striking and counter-intuitive example is provided by the Mpemba effect. Known already to Aristotle [1] but named after a high-school student E. B. Mpemba [2], the effect is commonly described as a curious phenomenon where initially prepared hot water freezes faster than cold water under otherwise identical macroscopic conditions, when both are cooled by the same cold bath. Due to the complexity of the phenomenon, the precise mechanism and conditions for the occurrence of the Mpemba effect have been under debate. Several explanations have been put forward to the particular mechanism for the Mpemba effect in water, highlighting possible roles for evaporation [3], supercooling [4], convection [5], particular properties of the hydrogen bonds [6,7], freezing-point depression by solutes [8] and a difference in the nucleation temperatures of ice nucleation sites between samples [9]. Moreover, the status of the Mpemba effect in water as an experimental finding has been recently called into question [10,11]. Indeed, subtleties of the liquid-solid transition make the precise definition of the effect difficult. For example, when does freezing occur? How well is the freezing point defined? Is a small left-over amount of vapor or liquid tolerated in the description? * mvucelja@virginia.edu; orenraz@gmail.com It is possible to view the effect as a particular example of a relaxation process far from equilibrium: the Mpemba effect is defined by a quenching process -cooling through a quick change in the ambient temperature, achieved by putting the system in contact with a new, colder, thermal bath. In contrast to quasi-static cooling, where the system is in equilibrium at each instant of the cooling process, quenching is, in general, a far-from-equilibrium process. Indeed, anomalous thermal relaxations are not unique to water, and similar effects have been observed in various other systems, e.g. magnetic alloys [12], carbon nanotube resonators [13], granular gases [14], clathrate hydrates [15], polymers [16] and even dilute atomic gas in an optical resonator [17].
Microscopically, the Mpemba effect occurs when the initially hotter system takes a non-equilibrium "shortcut" in the system's state space and thus approaches the new equilibrium faster than the initially colder system. A phenomenological description of such behavior was recently proposed by Lu and Raz within the framework of Markovian dynamics [18]. In this picture, a Mpemba like behavior can be studied in a large variety of systems (as many processes in physics and chemistry are Markovian [19]), and in particular, in small systems that can not be adequately described by macroscopic thermodynamics alone. An inverse Mpemba effect (associated with heating processes) can be described similarly. The suggested mechanism for the Mpemba effect raises several natural and important questions: (i) Does the mechanism require fine tuning of parameters, i.e., does it only occur in singular points of the model's parame-ter space? Is it robust to small changes in the system parameters? (ii) Does this mechanism survive the thermodynamic limit, or is it only a peculiarity of few-body systems such as those studied in [18]? One might intuitively expect that this mechanism does not apply to macroscopic systems, since in such systems the probability distribution is highly concentrated on specific points of the system phase space -those that minimize the free energy in equilibrium systems, and hence even if such "shortcuts" exist, the system cannot explore them.
The current manuscript contains the following contributions. First, using geometric insights on the relaxation dynamics in probability space we show that the Mpemba effect may be substantially enhanced on a discrete set of initial temperatures -a phenomenon we call the strong Mpemba effect. We show that these special initial temperatures can be classified by an integer, which we name the Mpemba index, and whose parity is a topological property of the system. Thus, the existence of a strong Mpemba effect is robust to small perturbations in the model parameters. Next, we study the effect in a thermodynamic system, focusing on a paradigmatic model: the meanfield Ising anti-ferromagnet, where a rich Mpemba-phase diagram is found. Using this model we demonstrate that even though in the thermodynamic limit the probability distribution is concentrated on specific points of phase space, the strong Mpemba effect still exists. Interestingly, we show that the strong Mpemba effect is tightly connected with another type of anomalous thermal relaxationthermal overshooting, in which the temperature of the system relaxes non-monotonically in time and overshoots the bath's temperature. We then provide an exact analytical calculation of strong Mpemba effect probability for an arbitrarily chosen set of energy levels in an idealized "isotropic" model. Comparison of these analytic results with the dynamics of the same set of energy levels with random barriers, gives a surprisingly good quantitative agreement. Lastly, we numerically study the strong Mpemba effect in the random energy model (REM) with random barriers, and find the scaling of the probability of the strong effect with the system size. This scaling suggests that again the effect can be observed in the thermodynamic limit.
The manuscript is organized as follows. In section II we give details on the explicit form of the Markovian dynamics, and in section III we define the strong Mpemba effect and describe its geometric meaning. Next, in section IV a study of the dynamics of an anti-ferromagnetic Ising model on a complete bipartite graph reveals a remarkable phase diagram exhibiting a variety of phases with various values of the Mpemba index, I M = 0, 1, 2 and phases with both direct and inverse strong Mpemba effect. In section IV D we show that the Mpemba effect also appears in the thermodynamic limit of the antiferromagnetic Ising model, and lastly, in section IV F we show that in this specific model the strong Mpemba effect implies overshoot in the temperature during relaxation. In section V A we study the probability of the occurrence of an odd Mpemba index for a system with random barriers taken from a very wide distribution and a fixed set of energies. For this purpose, we use as a model a statistically isotropic ensemble of second eigenvectors of the driving rate matrix. In particular, we find that the probability of a strong Mpemba effect is inversely proportional to the bath temperature, where the proportionality constant depends on the first few moments of the energy level distribution. In section V B we study the Mpemba effect in REM with random barriers.

II. SETUP AND DEFINITIONS
We consider Markovian relaxation dynamics, as given by the Master equation [20] where p = (p 1 , p 2 , ...), and p i (t) is the probability to be in state i at time t. Here, the off-diagonal matrix element R ij is the rate (probability per unit time) to jump from state j to i. The state i of the system is associated with an energy E i , and we focus on relaxation dynamics for which the steady state of Eq. (1)) is given by the Boltzmann distribution, where T b is the temperature of the bath, Z(T b ) = i e −β b Ei is the partition function at T b , and throughout the paper β ≡ (k B T ) −1 is the inverse temperature (in particular β b = (k B T b ) −1 ). Moreover, we also assume that the rate matrix R obeys detailed balance, and thus can be written in the form (see e.g. [21]): where B ij = B ji can be interpreted as the barrier between the states, and Γ is a constant with the proper units. At long times, the Markov matrix (4) drives an arbitrary initial distribution to the Boltzmann distribution associated with the bath temperature T b . Note that if the matrix R does not satisfy the detailed balance condition, its steady state does not represent an equilibrium since it has non-vanishing current cycles. Interestingly, a direct and an inverse Mpemba-like effects were recently discovered in driven granular gases where detailed balance is violated [14]. Although our approach may be useful also for such non-equilibrium steady states, for simplicity we limit our discussion to systems obeying detailed balance.
In the Mpemba effect scenario, the initial condition for Eq. (1) is the thermal equilibrium for some temperature During the relaxation process, the distribution p -i.e. the solution of Eq. (1) -can be written as where the rate matrix R has (right) eigenvectors v i and eigenvalues λ i , The largest eigenvalue of R, λ 1 = 0, is associated with the stationary (equilibrium) distribution π(T b ), whereas all the other eigenvalues have negative real part, 0 > Re λ 2 ≥ Re λ 3 ≥ ..., and they correspond to the relaxation rates of the system. The equilibration timescale is typically characterized by −(Re λ 2 ) −1 [22]. Any detailed balance matrix R can be brought to a symmetric formR via the similarity transformation, where F ij ≡ e β b Ej δ ij . The matrixR has the same eigenvalues as R, and has an orthogonal set of real eigenvectors. In particular This form will be useful in what follows.

The Mpemba effect
A simple criterion for the presence of a Mpemba effect for the relaxation process in Eq. (1) was given by Lu and Raz [18]. When | Re λ 2 | < | Re λ 3 | (namely when they are not equal), the probability distribution (6) can be approximated, after a long time, as In this case the Mpemba effect is characterized by the existence of three temperatures: hot, cold and the bath, The coefficient a 2 can be derived as follows: multiplying Eq. (6) with f 2 F 1/2 from the left, substituting v i = F −1/2 f i and using the fact that f i form an orthogonal basis, one get f 2 · F 1/2 p(T ; t) = a 2 (T )||f 2 || 2 e λ2t . Therefore for an evolution starting at a given initial probability p init we have that a 2 is the corresponding overlap coefficient between the initial probability and the second eigenvector f 2 : At the bath temperature this coefficient vanishes, a 2 (T b ) = 0 (as in this case F 1/2 p init = f 1 , which is orthogonal to f 2 ), and it increases in absolute value as the initial temperature departs from the bath. Therefore, to determine whether the Mpemba effect exists one has to look for non-monotonicity of a 2 (T ).
In the next section, we define the strong Mpemba effect, introduce an index to characterize the strong effect and describe the geometrical interpretation of the effect.

III. THE STRONG MPEMBA EFFECT, ITS INDEX AND ITS PARITY
Our first contribution is the observation that a stronger effect (even shorter relaxation time) can occur: a process where there exists a temperature T M = T b such that We call such a situation a strong direct Mpemba effect if T M > T b and a strong inverse Mpemba effect if T M < T b , as at T M the relaxation process is exponentially faster than for initial temperatures slightly below or above it.
Since there is essentially no difference between the direct and inverse effects, we refer to both of them as strong Mpemba effects. The strong Mpemba effect implies the existence of the "weak" effect, as in order to cross zero, a 2 has to be a non-monotonic function of temperature (because a 2 (T b ) = 0, whereas a 2 = 0 slightly above and below T b ) [24].
To study the strong Mpemba effect, we define the Mpemba indices as: and the total index as: I dir M changes its value when the number of zero crossing of the graph of a 2 (T ) changes in the interval T b < T < ∞. This implies that it is robust, as discussed in section III B.

A. The geometry of the strong Mpemba effect
The geometry of the problem is schematically illustrated in Fig. 1 for a three state system. The set of all points p = (p 1 , p 2 , p 3 ) which are normalized (p 1 + p 2 + p 3 = 1) and non-negative (p i ≥ 0) form the probability simplex -the blue triangle in the figure. The set of all Boltzmann distributions π(T ) form the Boltzmann curve -the red line, which has two boundaries: the π(T = 0) where the probability is concentrated at the lowest energy state (blue point), and π(T = ∞) -the maximally mixed state in the middle of the simplex where all the states are equally probable (the red point). The set of all p for which a 2 = 0 is illustrated by the intersection of the green plane with the blue triangle.The Boltzmann curve intersect the a 2 = 0 plane at π(T b ) (pink point) since, being the equilibrium distribution at T b , a 2 (T b ) vanishes by definition. However, in the specific example, the two boundaries of the Boltzmann curve are both on the same side of the a 2 = 0 hyperplane, therefore there must be another point -marked by the green point in the figure -at which the Boltzmann line cross the a 2 = 0 hyperplane again. This point corresponds to T M where there is a strong Mpemba effect. Topologically, having co-dimension 1, the a 2 = 0 hyperplane separates the probability simplex into two disjoint sets. The parity of the number of times a continuous curve crosses this hyperplane depends only on the two boundaries of the curve: if they are both in the same set, then the number of crossing is even -and hence the parity of I M is odd (as it does not count the crossing at T b ), and if they are in a different set then the number of crossing must be odd, with an even parity for I M .
FIG. 1. The geometry of the strong Mpemba effect in a threestate system. The probability simplex is illustrated by the blue triangle. The set of all equilibrium distributions form the red curve. The blue and red points on the curve correspond to T = 0 and T = ∞. The Mpemba index is non-zero if the equilibrium curve crosses the a2 = 0 plane (green) not only at the bath temperature (illustrated here by the pink point) but also at some other temperature (the green point).
Given T b , a sufficient condition for the strong direct Mpemba effect to occur is obtained by determining whether a 2 changes sign going from a 2 (T b + ε) to a 2 (T = ∞). This can be expressed by where we used a 2 (T b +ε) ≈ ∂ T a 2 | T =T b ε and θ is the Heaviside step function. The argument of the step function in Eq. (15) is positive if a 2 (T ) has an odd number of zero crossings, thus Eq. (15) described the parity of the number of zeros. In particular, if P(I dir M ) = 0, we are assured to have at least one crossing, and so P(I dir M ) serves as a lower bound on the number of initial temperatures for which the direct strong Mpemba effect occurs. Similarly, the parity for the strong inverse Mpemba effect is and the parity of the strong Mpemba effect is As already mentioned above, in some situations there are zeros of a 2 (T ) that are not accompanied with a sign change. This happens when a 2 (T M ) = 0 and a 2 (T M ) = 0 simultaneously. Such points appear on the boundary between areas in parameter space with I M = 0 and I M = 2, e.g. on the line separating purple and green areas in Fig. 2. In these cases, P(I M ) is no longer the exact parity, but still serves as a lower bound to the number of crossings. Fig. 1 provides a helpful three-dimensional picture for the strong Mpemba effect in a three-state system. In fact, the existence of the strong Mpemba effect is always essentially a three-dimensional problem, when projected to the proper sub-plane: as discussed above, the strong effect can be deduced from the relative directions of the following three vectors: (i) The tangent to the equilibrium line at the bath temperature; (ii) The vector connecting π(T b ) and π(T = ∞) (or π(T = 0) for the inverse effect); and (iii) the direction of the slowest dynamic, v 2 defined in Eq. (7)). This observation plays a crucial role in Sec.(V A), where we estimate the probability to observe the strong Mpemba effect in a class of random models.

B. Robustness of the Mpemba Index
The above geometric interpretation implies that the Mpemba index is a robust quantity, as we discuss next. Consider a small perturbation of order in the physical quantities, i.e. in the energies, the barriers and the temperature of the bath. The corresponding change in R is also of order . This perturbation in R changes both p init (T ) and f 2 , and hence by Eq. (11), also the graph of a 2 (T ) in the relevant interval of T . But it does not change π(T = ∞), which is always the maximally symmetric point. Similarly, π(T = 0) is the lowest energy point, which changes only if the perturbation changes the ground state by an energy level crossing.
For the perturbation in the physical parameters to changes the number of zero crossings, one of the following cases has to occur: (i) p init can change abruptly even with a small perturbation in R, if there is a first order phase transition in the system and therefore the equilibrium distribution changes discontinuously. An example for such a case is discussed in Sec. IV D. Note that for this to happens, the perturbation in the relevant parameters has to be large enough compared to the distance from the value at which there is a first order phase transition. (ii) f 2 can change abruptly when the perturbation changes the order of λ 2 and λ 3 , namely the order of the eigenvalues changes and therefore the direction of the eigenvector jumps. As the perturbations of the eigenvalues are of order too, the spectral gap λ 2 − λ 3 defines the stability region in which these changes are not expected. In other words, if the perturbation in R is small compared to λ 2 − λ 3 then such a jump is not expected. (iii) A small perturbation in R can cause two zeros to "annihilate" each other in a saddle-node bifurcation, and similarly two new zeros can be generated. However, in these cases the parity of I M does not change. (iv) A zero can move through T b (this is analogues to a transcritical bifurcation). In such a case both I dir M and I inv M change by one, but the parity of I M does not change. Lastly, (v) a zero can "vanish" in the T → ∞ limit or at T = 0. But as discussed above, the end points π(T = 0) and π(T = ∞) do not move, so this can only happen if f 2 changes its direction.
From the five cases discussed above we can conclude that the parity of I M can change only if there is a phase transition in the system, λ 3 becomes larger than λ 2 or if the perturbation changes the sign of one of a 2 (0) and a 2 (∞). Therefore, the parity is stable in some range, which depends on details as the spectral gap λ 2 − λ 3 , the distance of parameters from a phase transition and the angle between π(T = 0), π(T = ∞) and f 2 .
The above argument for stability can be explained using the geometric picture for the Mpemba effect discussed above. In any system, the a 2 = 0 hyperplane separates the probability space into two disjoint sets as discussed in the three state system above. The Boltzmann curve intersects the a 2 = 0 hyperplane at T = T b . Any additional temperature T M for which a 2 (T M ) = 0 is an intersection between the curve and the hyperplane, and I M counts these additional intersections. Topologically, the parity of the number of crossing between a continuous curve and a hyperplane of codimension 1 depends only on the boundaries of the equilibrium curve, namely in our case on π(T = 0) and π(T = ∞). If they are both on the same side of the hyperplane then the number of crossing is even, and if they are on different sides then the number of crossing is odd.
To appreciate the above topological aspect, let us contrast it with a (possibly only hypothetical) "super strong Mpemba effect", where there exist a temperature T SM at which a 2 (T SM ) = a 3 (T SM ) = a 4 (T SM ) = 0. In other words, consider a case where the coefficients of π(T SM ) vanish along both v 2 , v 3 and v 4 . This implies an even faster relaxation than the strong effect. The condition for this super strong effect, a 2 (T SM ) = a 3 (T SM ) = a 4 (T SM ) = 0, defines a hyperplane of co-dimension 3, that can intersect the probability simplex (which is of co-dimension 1) in a co-dimension 2 hyperplane. This hyperplane does not separate the space into two disjoint sets, and the topological argument does not work anymore. An example for this fact, consider a codimension 2 hyperplane in a 3D space, which is a straight line. It does not separate the 3D space into two disjoint sets as the a plane does. Now consider an equilibrium locus that crosses the super-strong straight line in a 3D probability simplex. A small perturbation in the model parameters deforms a bit the equilibrium locus, and generically separates the equilibrium locus from the straight line. Therefore, even an infinitesimal perturbation can (and usually does) destroy the super-strong Mpemba effect. In this sense, the strong Mpemba effect has a topological stability, but the super strong effect does not. We therefore do not expect to observe the super strong effect in systems which are not fine tuned.

IV. MEAN FIELD ISING ANTIFERROMAGNETIC MODEL WITH GLAUBER DYNAMICS
The mechanism for the Mpemba effect suggested in [18] was so far demonstrated only in microscopic systems with a few degrees of freedom. However, all the experimental observations of similar effects are in macroscopic systems, with a huge number of micro-states. To discuss the applicability of the mechanism in such macroscopic systems, we next consider the Ising model, with anti-ferromagnetic interactions and mean-field connectivity, on a complete bipartite graph. This is a classical many-body model which has been studied extensively, and whose phase diagram can be calculated exactly (see figure 6 in [25]). As described below, this model shows a rich Mpemba behavior, which survives the thermodynamic limit.

A. The Model
In mean-field models, each spin interacts equally with all the other spins in the system. To generate a model of antiferromagnetic interactions in the mean field approximation, we consider a system with a total number of N spins, half of them on each "sub-lattice" or subgraph. Each spin interacts antiferromagnetically with all the spins in the other sub-graph, but spins on the same sub-graph do not interact at all. The interaction strength between the spins is fixed. This type of interaction can lead to an "antiferromagnetic phase" in which the spins in one sub-lattice are predominantly in the up state, while most spins in the other sub-graph point down.
Let N 1,↑ , N 1,↓ (N 2,↑ , N 2,↓ ) be the number of spins pointing up and down on sub-graph 1 (sub-graph 2). We define the two magnetization densities on sub-graphs 1 and 2 as: Although the system has 2 N different microstates, all microstates that correspond to the same values of N 1,↑ and N 2,↑ are equivalent, since the interaction strength is "position" independent (mean field). Thus, the Hamiltonian of this model is only a function of x 1 and x 2 and is given by where J is the coupling constant, H is the external magnetic field and µ is the magnetic moment. In the antiferromagnetic case, the coupling constant is negative, J < 0, and for simplicity we choose the units such that J = −1 and µ = 1 . The dynamics we consider for this model is Glauber dynamics, with only a single spin flip transitions allowed. Under this assumption, the rates of flipping a spin up or down in sub-graphs 1 and 2 are given by where R u1 (x 1 , x 2 ) is the rate of flipping a spin up in subgraph 1 and R d2 (x 1 , x 2 ) is the rate of flipping a spin down in sub-graph 2. The numerators in Eqs. (20) are the combinatorial factors that take into account how many spins can be flipped in the specific state of the system, and the denominator is the standard Glauber factor, 1/(1 + e β b ∆E ), where ∆E is the difference of energies before and after the spin flip [26].

B. Mpemba Index phase diagram
The Mpemba-index phase diagram of this model was calculated numerically for N = 400, and is shown in the upper panel of Fig. 2. At each point in the figure, that is -for each temperature T b and magnetic field H of the environment, we calculated (numerically) the coefficient a 2 (T ) of the slowest relevant eigenvector of the corresponding Glauber dynamics (Eq. 20) at each point along the equilibrium line. From the monotonicity and zero crossing of these coefficients, a 2 (T ), we have deduced what types of Mpemba effect exist at this point. The phase diagram in Fig.(2) is quite rich and has 8 different phases, differentiated through their colors, including regions with odd and even Mpemba index, existing for the direct inverse or both effects.
To make sure that the observed phase diagram is not dominated by the number of spins in the system, we repeated this calculation with N = 70 and checked that the phase diagram looks essentially the same (Fig. 2, lower  panel). Moreover, we have checked other form of rates -Metropolis and heat-bath dynamics, both with single spin flips only. Although the exact locations of the different phases are not identical in the different dynamics, the main features in the phase diagram are similar in all of them. An example for such a feature is the line at H = 1 across which the Mpemba phase changes. To explain this feature, we next consider the thermodynamic limit of this model.  [25]. Lower panel -the same calculation for N = 70 spins. Although the exact location of boundaries between the different phases is not identical in the two computations, the overall structure is the same. The three dashed lines are lines of constant magnetic field H = 1.03, H = 1 and H = 0.97, and their corresponding equilibrium loci in the thermodynamic limit are given in Fig. 3. Note that the phase diagram has an abrupt jump at H = 1, which corresponds to the jump in the equilibrium locus in Fig. 3.

C. The Thermodynamic Limit
Let us take the thermodynamic (N → ∞) limit for the mean field anti-ferromagnet model described above. To this end, we first write explicitly the master equation using all the single flip rates. At the configuration (x 1 , x 2 ), a single spin in each sub graph can either flip from up to down or from down to up. Therefore, there are four different terms in the master equation corresponding to leaving the current configuration, and similarly four transitions into the specific configuration: where ∆x is the change in the variable x due to a single spin flip. In the limit N → ∞ we approximate x 1 and x 2 as continuous variables. Expanding both p and all the terms of R to first order in ∆x we get a Fokker-Planck like equation Note that in this case there is no diffusion, as the corresponding term vanishes in the N → ∞ limit. Hence, it originates from a Langevin equation without random noise, namely from a deterministic equation for x 1 and x 2 . For such deterministic motion, an initial distribution which is a delta function stays a delta function at all times, and it is therefore enough to know the evolution of the averages Using these definitions, we write an "equation of motion" for the averages of x 1 and x 2 by substituting the values of the rates in Eq. (20) into Eqs. (22,23). After some algebra these give: Unfortunately, these equations are not always linearly stable: for some values of H, T b , x 1 , x 2 a small perturbation in the initial values of x 1 , x 2 changes significantly the trajectory. For example, when the initial condition has x 1 = x 2 the symmetry of the dynamic keeps x 1 and x 2 equal at all times, even if the equilibrium distribution which corresponds to the specific T b and H is different. In such cases, any infinitesimal perturbation (in the initial condition or during the dynamic) results in relaxation towards the equilibrium, rather than following the solution of the above equations. Fortunately, in a large fraction of the parameter space (H, T b ), as well as in the vicinity of all the fixed points, the above equations are stable. When stable, these non-linear equations describe the temporal evolution of the macroscopic system, and we can use them to understand the Mpemba behavior of the system. Using the above result, let us look at the equilibrium locus in the thermodynamic limit. For each value of H and T b , the equilibrium values of x 1 , x 2 , which we denote by ξ 1 , ξ 2 , are the steady state Eq. (25), namely the solution of For each value of H, the equilibrium line can therefore be found using Eqs. (26) to calculate ξ 1 (T b ) and ξ 2 (T b ), for 0 ≤ T b ≤ ∞. Note that Eqs. (26) are symmetric to exchanging x 1 and x 2 . We therefore limit ourselves, without loss of generality, to x 1 ≤ x 2 . Examples for the equilibrium locus for H = 0.99, H = 1 and H = 1.01 are shown in Fig. 3, where for each T b we have numerically found the equilibrium by solving Eq. (26). As expected, at T b → ∞ both ξ 1 → 0 and ξ 2 → 0. Importantly, the equilibrium line is not a simple convex line, therefore increasing the distance along the line does not necessarily increase the changes in the magnetizations. Moreover, at H = 1 the equilibrium locus has a singular transition demonstrated in Fig. 3. The sharp change in the equilibrium line at H = 1 corresponds to the first order phase transition in the model at H = 1 when T b = 0. This transition can be demonstrated by considering the limit T b → 0: for H > 1, the arguments of the hyperbolic tangents in Eq. (26) approach +∞ asymptotically in the limit T b → 0, and hence ξ 1,2 → 1. In contrast, for H < 1, ξ 1 → 1 and ξ 2 → −1 are the asymptotic solutions at T b → 0. As discussed in Sec.(III), this sharp transition in the shape of the equilibrium line naturally corresponds to the sharp transition in the Mpemba-phases in Fig. 2: when the equilibrium locus abruptly changes, so does the coefficient along the slowest relaxation mode, a 2 (T ). Although this argument in principle should work only in the thermodynamic limit, it is evident from Fig. 2 that in practice it works well already at N = 70.
So far we have seen that some features of the finite system Mpemba phase diagram can be explained using the thermodynamic limit. In the next section we discuss the existence of Mpemba effects in the thermodynamic limit, and their relations to the finite N system.

D. Weak and Strong Mpemba effects in the Thermodynamic limit
For systems with a finite number of states and a given set of environmental parameters, we have a simple prescription to check what types of Mpemba effects exist: the monotonicity (weak effect) and zero crossings (strong effect) of the coefficient along the slowest dynamic, a 2 (T ), encapsulate all this information. In the thermodynamic limit, we cannot use the same method, as rarely does the coefficient a 2 (T ) have an analytic expression at finite N , for which we can take the thermodynamic limit, N → ∞. Likewise, the direct calculation of a 2 (T ) in the infinite system is often not feasible. This is somewhat unfortunate, because all the experimental observations of Mpemba effects mentioned above are in macroscopic systems, hence it is not clear that the mechanism we suggest is relevant for such systems. Moreover, the existence of the effect in macroscopic systems does not follow trivially from its existence in small systems. Although phase space becomes larger in large systems thus more shortcuts may exist, in the thermodynamic limit the probability distribution is concentrated in a tiny portion of the systems phase space, which suggests that the system would rarely explore the extended phase space, and therefore these shortcuts might not be as relevant.
Although we cannot use a 2 (T ) to analyze the existence of Mpemba effects in the thermodynamic limit of the antiferromagnet Ising model, for any environmental conditions (T b , H) and two temperatures T h and T c , Eq. (25) can be used to compare the relaxation trajectories ini-tiated from the corresponding equilibrium distributions. A natural and physically motivated distance function in this case is the free energy difference between the current state and the equilibrium state, namely (27) where the free energy is given by [25]: and H(x 1 , x 2 ) is given by Eq. (19)). If the initial condition with a longer distance from the equilibrium becomes, after some finite time, closer to equilibrium, then we know that there is a Mpemba effect in this system. However, checking if a Mpemba effect exists at a point using this approach is tedious -it requires solving the relaxation trajectories for all initial conditions.Luckily, checking if strong effects exist and identifying their index is a much easier task. To this end we can linearize Eqs. (25) near the equilibrium point corresponding to (T b , H). Denoting the differences from the equilibrium by ∆x i = ξ i − x i , we can write for small ∆x i These linearized equations have two relaxation eigendirections -a fast direction and a slow direction, with relaxation rates given bȳ Unless the initial condition is such that at the long time limit the coefficient of points on its trajectory along the slow direction (namely the eigenvector corresponding tō λ 1 ) is zero, at long enough time the relaxation is from the direction corresponding to the slow direction. The number of trajectories that start on the equilibrium locus and approach the equilibrium point asymptotically from the fast direction is the Mpemba index, and the corresponding initial conditions show strong Mpemba effect. To find if such initial conditions exist, we can shoot backwards in time a solutions to Eqs.(25) that approach the equilibrium from the fast direction (there are two such trajectories -one from each side of the equilibrium locus). The number of crossing between these shoot-back trajectory and the equilibrium locus is the Mpemba index.
As an example, consider the relaxation dynamic for environment with H = 1.1 and T b = 0.5. The equilibrium locus as well as the relaxation trajectories from different initial temperatures are plotted in the upper panel of Fig. 4. As can be seen in the figure, there is a "fast" and a"slow" direction to the relaxation process, and essentially all the trajectories relax to equilibrium from the "slow" direction -except for a single trajectory (the red dashed trajectory) that relax directly from the fast direction. This special trajectory is the strong inverse Mpemba with I M = 1 -in agreement with the finite state phase diagram in Fig. 2 for these values of H and T b . In this case, the ratio between the relaxation ratesλ 1 andλ 2 given in Eq. (31) is 14.7. In other words, not only the strong Mpemba initial condition relaxes exponentially faster, the relaxation rate is an order of magnitude higher that from any other initial temperature. Indeed, as shown in Fig. 5, the strong effect trajectory relax significantly faster towards the equilibrium compared to all the other initial temperatures. E. Comparing the thermodynamic limit with a finite N system size So far we have seen the strong Mpemba effect in both the finite state N = 400 and the thermodynamic limit of the anti-ferromagnetic mean field Ising model. However, it is not obvious that the two effects are trivially related: In the thermodynamic limit the effect was derived by linearizing the non-linear equations for the order parameters, Eqs. (25), not by considering a 2 (T ) in the N → ∞ limit, which is intractable. One way to compare the two cases is to calculate a Mpemba phase diagram for the thermodynamic limit, and compare it with Fig. 2. However, in the paramagnetic phase (x 1 = x 2 ) we cannot use Eqs.(25) as they are not linearly stable, as discussed above. Therefore, we performed instead several other comparisons as presented next.
One hint that the two effects are nevertheless related is given by the sharp transition of the equilibrium line at H = 1 and the corresponding jump in the Mpemba index shown in Fig. 2 and discussed above. To further convince that the strong Mpemba in the thermodynamic limit corresponds to the finite system case, the temperature at which a strong effect occur, T M (N ), was calculated for various values of system size N . Fig. 6 shows that indeed T M (N ) converges, at the large N limit, to the temperature T ∞ M at which a strong effect occur, in the thermodynamic limit, for the chosen T b and H. Similar behavior was observed in a wide range of temperatures and magnetic fields.
Additional comparison between a finite system and the thermodynamic limit is shown in the lower panel of Fig. 4. For N = 400, we have calculated the equilibrium distribution as a function of the temperature, and "projected" it into the (x 1 , x 2 ) plane by calculating the equilibrium averaged magnetization in each sublattice. This   (20). These trajectories in probability space were projected to the (x 1 , x 2 ) plane, and are shown by the colored thin lines. Lastly, the trajectory of the strong effect T M = 0.148 that solves a 2 (T M ) = 0, was calculated too, and is plotted by the thick blue dashed line. For comparison, the strong Mpemba trajectory in the thermodynamic limit, initiated at T ∞ M = 0.1377, is plotted in thick dashed red line. As shown in the inset, all the trajectories except the strong Mpemba approach the equilibrium point (the blue circle) from the "slow" direction, and the strong effect approaches the equilibrium point from a different -"fast" direction. Note that the equilibrium line and the trajectories in the thermodynamic limit are not identical to those of the finite system, hinting that N = 400 is not large enough to match the thermodynamic limit. Nevertheless, this example demonstrates that the a 2 (T M ) = 0 in finite size systems maps into the strong Mpemba mechanism in the thermodynamic limit. Although not a mathematical proof, this analysis hints that the strong Mpemba effect in the thermodynamic limit is a consequence of the strong effect in the finite system.  6. Comparison between the strong effect in N spin system and the strong effect in the thermodynamic limit. The temperature from which there is a strong effect at finite N system, TM (N ), converges to the temperature from which there is a strong effect in the thermodynamic limit, T ∞ M .

F. Temperature overshooting during Relaxation
As discussed above, the existence of a Mpemba effect can be checked by calculating the distance from equilibrium as a function of time from two different initial points on the equilibrium lines. This requires a choice of some reasonable distance function, e.g. the free energy distance (see Eq. (27)). In the specific case of the mean field anti-ferromagnet model there is another natural option: even though the system is not in equilibrium through the relaxation, it is possible to associate a temperature with each state during the relaxation process, and use this temperature to compare different relaxations. This temperature does not have all the properties commonly required from a distance function, e.g. it is not monotonically decreasing in a relaxation, nevertheless, it shed light on additional counter-intuitive aspect of thermal relaxations far from equilibrium: as shown below, the temperature can overshoot the environment temperature. In other words, a hot system, when coupled to a cold bath, can reach during its relaxation process temperatures which are lower than the environment's temperature.
To associate a temperature for each point in the relaxation process, let us use the following coordinate transformation, from (x 1 , x 2 ) to (T eq , H eq ), defined by: The physical significance of this transformation can be understood by a simple algebraic manipulation of the above equations that gives Eq. (26). Comparing these to Eq. (25) one notes that for environment with temperature T b = T eq and external magnetic field H = H eq , the specific state given by (x 1 , x 2 ) is the equilibrium. In other words, if during the relaxation process, when the system is in the state (x 1 , x 2 ), the system is decoupled from the current environment and coupled to a different environment with T b = T eq and H = H eq , then the system would be in equilibrium with the new environment. It is therefore natural to interpret H eq and T eq as the magnetic field and temperature of the system itself.
Before proceeding, two comments on the above mapping are in place. (i) Note that the transformation is singular at x 1 = x 2 as the denominator in Eqs. (33) vanishes. In other words, we cannot associate a single temperature and magnetic field for states in which x 1 = x 2 .
(ii) The ability to associate equilibrium temperature and magnetic field to most states of the system is a very nongeneric property. It is a consequence of the fact that the number of parameters in the model is identical to the number of order parameters describing the system in the thermodynamic limit. Luckily, in the thermodynamic limit of this model the probability distribution becomes a delta function with exactly two order parameters.
Using the mapping in Eq. (33), we plotted in Fig. 7 the temperature of the system as a function of time for various initial conditions along the equilibrium line. As can be seen, not only the temperature curves cross -namely a Mpemba effect occurs, but also for some relaxation trajectories the temperature is non-monotonic as a function of time. Moreover, systems that were initiated at temperatures lower than the environment's temperature can reach temperatures which are higher than that of the environment in their relaxation. Similar non-monotonic relaxations were discussed in the context of non-Markovian thermal relaxation [27] or finite baths [28], but as far as we know this is the first example for such non-monotonic relaxation in Markovian dynamics and in the thermodynamic limit.
It is interesting to note that the temperature overshoot is tightly connected with the strong Mpemba effect. To explain this, let us examine Fig. 4 carefully. Initial temperatures to the left of the strong Mpemba relaxation trajectory approach equilibrium from one side of the slow direction, and initial temperatures to the right of the strong Mpemba trajectory approach the equilibrium from the opposite direction, which is also a slow direction. Opposite directions mean opposite directions in the coordinate T eq , namely there are trajectories that approach the equilibrium both from higher and lower temperatures compared to the environment. Conversely, if there are trajectories that approach equilibrium from both higher and lower temperature, they must approach equilibrium from opposite directions of the slow relaxation. Therefore, by continuity there must also be a trajectory that approach equilibrium from the fast direction -and this trajectory corresponds to a strong effect. FIG. 7. Teq − T bath as a function of time for various initial temperature along the equilibrium line. The inverse Mpemba effect is shown here as a crossing of two curves -initially colder system heats up faster. The figure also shows that the temperature can overshoot the environment temperature -the system can reach equilibrium temperatures which are higher than that of the environment.

V. HOW GENERIC IS THE MPEMBA EFFECT IN THE REM MODEL?
So far, we considered the Mpemba effect in a specific model -the mean field antiferromagnet Ising model. Is the existence of Mpemba effect a special property of this model, or should we expect to see similar effects in many other models? To address this question, we next evaluate the probability to have a Mpemba effect in a class of models with random parameters. As discussed below, the strong Mpemba effect plays a crucial role in our ability to estimate the probability of having a Mpemba effect.

A. Analytical estimates -The Isotropic ensemble
A proper analysis of the probability to find a Mpemba effect in classes of random models is a formidable challenge: one must perform the rather difficult calculation of the second eigenvector f 2 and the coefficient a 2 in Eq. (11) as a function of the initial temperature, the energies and the barriers and then average over the ensemble. Even the simpler problem of analyzing the strong Mpemba effect requires facing the daunting task of analyzing the number of zeros in Eq. (11).
To gain some analytical insight, we proceed by estimating the strong Mpemba probability in an ensemble of relaxation dynamics which we call the isotropic ensemble. The ensemble is chosen to represent a wide distribution of barriers, so that the distribution of eigenvectors of the relaxation modes is as isotropic as possible consistent with a given target thermal distribution. Explicitly, given a set of L energies {E 1 , E 2 , ..., E L }, we average over an ensemble of random f 2 eigenvectors that are orthogonal, in the sense of the quadratic form in Eq. (11), to the equilibrium distribution with a given bath temperature T b .This approach allows us to perform analytically the ensemble averaging. We then compare our analytic results for the isotropic ensemble with direct numerical calculations on the matrix Eq. (4) with fixed energies and random barriers, and find a surprisingly good agreement in certain parameter regimes.
For a given set of energies E 1 , ..., E L and dynamics prescribed by the Markov matrix Eq. (4) the steady state distribution is the Boltzmann distribution at T b : π(T b ). The first eigenvector of the symmetrized Markov matrix R, is andRf 1 = 0. The second eigenvector ofR, f 2 , together with the initial condition π(T ), determine the coefficient a 2 , which according Eq. (11) is We obtain the explicit expression for the parity of the direct Mpemba index as a function of f 2 by plugging Eq. (35) into Eq. (15) and find where is the average energy in equilibrium at T b . We can represent Eq. (36) in the form with the vectors u dir and w defined as Note that the vectors u dir and w appearing in this form depend solely on the set of energies and on the bath temperature -they are independent of the barriers. Moreover the form of Eq. (37) has a simple geometric meaning. To see it we single out the components of f 2 in the plane spanned by the (non-orthogonal) vectors u dir , w.
Choosing f as the component of f 2 parallel to u dir , we have: + terms orthogonal to u dir and w. (40) In this basis (f 2 · u dir )(f 2 · w) is equal to Therefore on the f , f ⊥ plane, the region satisfying P(I dir M ) = 0 is a double wedge (see Fig. 8). The boundary of the region is associated with the lines f ⊥ = −f /K and f = 0. The same treatment is also possible for the inverse Mpemba effect. For example, assuming, for simplicity, a non-degenerate ground state, and ordering the energies so that E 1 is the ground state energy, we find that P(I inv M ) is given by Eq. (37) with the replacement Next, we formulate the averaging over the admissible f 2 vectors. In the class of random relaxations we consider, we generate f 2 by picking a random vector g = (g 1 , ..., g L ), and obtaining from it a random vector orthogonal to f 1 (by subtracting the projection of g on f 1 ) The distribution of the g vectors is taken to be isotropic and therefore the projection of the distribution of gs onto the hyperplane orthogonal to f 1 is also isotropic. For this purpose, we take the g i s in g to be IID Gaussian variables. Analogously to the derivation of, Eq. (37), we plug Eq. (45) into Eq. (36) and separate the g i components.
We find that the direct Mpemba parity, for a particular realization of g i , can be written as with w is defined in Eq. (39) and u dir iso given by where L is the number of energy levels (or the system size). As before, we break g into the component parallel and perpendicular to u dir iso and find, as before: where K is defined in Eq. (42). The Gaussian IIDs g i have a rotationally-invariant joint distribution function, and therefore in any coordinate system the components are corresponding Gaussian IIDs. It follows that g , g ⊥ are Gaussian IIDs and have a rotationally-invariant distribution on the g , g ⊥ plane. On this plane, the region satisfying P(I dir M ) > 0 is a double wedge (c.f. Fig. 8), and the probability of g , g ⊥ to fall inside the wedge only depends on the wedge angle.
Geometrically, if φ is the angle between u and w, then Prob(P(I dir M ) > 0) = φ π when (u · w) > 0 (and Prob(P(I dir M ) > 0) = 1 − φ π when (u · w) < 0). Expressed explicitly in terms of u dir iso , w we find: To recap, the formula Eq. (50) represents, for a given set of energies {E i } and bath temperature T b , the probability that the direct Mpemba index is odd. Eq. (50) can be simplified for hot bath temperatures, k B T b max({E 1 , ..., E L }), and asymptotically it gives Here the constant C E depends only on the first few moments of the energy level distribution (for the explicit expression see the appendix VII). Fig. 9 shows a comparison of Eq. (50) with a random realization of L = 15 energies, and Mpemba index averaged over 4000 realizations of random barriers. The expression seems to capture surprisingly nicely the behavior of for a random draw of energy levels when the barriers distribution is wider than the distribution of energies, and the temperature is higher than the characteristic energy spread.
Similarly, for inverse Mpemba we have with w is defined in Eq. (39) and u inv iso given by where we assumed that E 1 is the lowest energy. As before, Substituting for w and u inv iso from Eqs. (39) and (53) we get The above expression simplifies in the limit of a very low bath temperature, k B T b (E 2 − E 1 ). Without loss of generality we set E 1 = 0 and obtain .
Simplifying the expression even further we get which is expected as there is no Mpemba effect for a two level system. It is important to note that the isotropic ensemble, while introduced for the purpose of enabling analytical averaging, is consistent with the assumptions of our relaxation dynamics. Namely, one can prove the following theorem.
Theorem: Given any choice of a real vector f 2 orthogonal to f 1 in (34), there exists a set of barriers B ij with relaxation dynamics obeying detailed balance (4) having F −1/2 f 2 as its slowest relaxation eigenvector.
The proof can be found in the appendix VIII. In the previous section we analyzed an analytically tractable model, however, in general, estimating the probability of the Mpemba effect is a daunting and often impossible task. In what follows, we numerically study the probability of having a strong Mpemba effect in the random energy model (REM). The random energy model was introduced by Derrida as an extreme limit of spin glasses [29]. It is the simplest model of a system with quenched disorder that has a phase transition. In the REM, L energy levels are IID random variables. The conventional choice for the probability distribution of E j 's is a Gaussian distribution where in order to have extensive thermodynamic potentials the variance depends on the system size. At temperatures lower than T critical ≡ σ E /(k B √ 2 ln 2) the system is trapped in a few low-lying states; this condensation phenomena is a phase transition and at the transition the free energy is non-analytic.
Note that quenched disorder is not necessary for the Mpemba effect. Our previous examples, such as the mean field Ising antiferromagnet and other examples discussed in [18] are a proof that quenched disorder is not an essential feature for this effect.
The Mpemba effect is a property of the system and its dynamics, thus to study it we need to specify the barriers B ij in (4). Here we chose B ij as IID random variables obeying a "truncated" Gaussian distribution and θ is the Heaviside step function. This particular choice of barriers can only impede the transition rates R ij , as in this case e −β b Bij < 1. The variance of the barriers is scaled with the system size like that of the energies, so that their ratio is system size independent. Note that numerous other choices of the dynamics for the REM have been studied in the past, most notably single spin flip dynamics (see e.g. [30][31][32] and references therein) and it would be interesting explore for Mpemba effects those other choices of REM dynamics as well.
Numerically we studied the parity of the direct Mpemba effect (see Eq. (15)), by exact diagonalization of an ensemble of REM with random barriers R matrices. As an example of typical numerical results, see Fig. 10, where L = 10 energy levels were chosen from a Gaussian distribution Eq. (59) and barriers chosen from Eq. (60). The bath temperature in the numerics was k B T b = 0.1 and k B T b = 1.0. Each data point was averaged over 10 5 realizations. From the ample numerical evidence we infer that the Mpemba effect occurs with finite probability, especially for T b < T critical case (left panel Fig. 10).
We also studied the system size dependence of REM with random barriers, see Fig. 11. The system size was varied (L ∈ [4,20]) and we took the bath temperature to be k B T b = 0.1. The energies were IIDs from Eq. (59) with variance σ 2 E log 2 L where σ E = 1.0 and barriers were IIDs from Eq. (60) with variance σ 2 B log 2 L. Each point on the density plot was averaged over 2 × 10 5 realizations. We notice that the probability of the parity being positive for the direct Mpemba index seems to be converging to a limiting value with increasing system size. Although, we tested small sizes, the convergence suggests the thermodynamic limit behavior.

VI. DISCUSSION
The Mpemba effect is a "short-cut" in relaxation time. The direct Mpemba effect implies that initiating the system at a particular hot temperature results in cooling down which is faster than any colder temperature, when the system is coupled to a cold bath. Possibly even more counter-intuitive is the inverse Mpemba effect where an analog effect happens in heating. Similarly to the direct Mpemba effect, in annealing one first heat the system and then cools it in a controlled manner such that it acquires desirable features (relaxes to the ground state, has fewer defects, etc.). More specifically, simulated annealing is a probabilistic technique used to find ground states [33][34][35], while annealing in metallurgy is used to make materials with larger mono-crystal domains and fewer defects [36]. It would be interesting to explore the connection between annealing and the Mpemba effect.
Markov Chain Monte Carlo (MCMC) algorithms are essential numerical tools broadly used in many branches of science to estimate steady-state properties of various systems [37]. It is often desirable to speed up the relaxation of a MCMC to the steady state, see, e.g., [38,39]. Our results serve as a proof of principle that in specific systems one could devise additional transition barriers (B ij s) that would cause speed up of a MCMC algorithm's relaxation to equilibrium by creating a strong Mpemba effect.
Approach to equilibrium often has a non-trivial relationship with the energy landscape and nature of the barriers. This is especially true in glassy materials and complex many-body systems. The approach to equilibration can even be used to explore structures in glassy systems and many-body systems experimentally [40]. One of the future directions is to deepen the understanding of the relation of the Mpemba effect and to the plethora of Lower bound for the probability of the strong Mpemba effect -the probability of the parity being positive for the direct Mpemba index (see Eq. (15)) for the case of REM with random barriers and different system sizes L ∈ [4,20]. The bath temperature is kBT b = 0.1. The energies were IIDs from Eq. (59) with variance σ 2 E log 2 L and σE = 1.0 and barriers were IIDs from Eq. (60) with variance σ 2 B log 2 L. Each point on the density plot was averaged over 2 × 10 5 realizations. We notice that the probability of the parity being positive for the direct Mpemba index seems to be converging to a finite limiting value with increasing system size.
nontrivial cooling phenomena present in glassy materials; such as memory, aging, and rejuvenation.
where K is given in Eq. (42). By plugging in Eqs. (38) and (39) in Eq. (42) we get At the high temperature limit, T b → ∞, we can expand K in small β b . To get the correct result, we have to expand all terms in the argument for the square root up to order β 2 b . Using arctan 1 K ∼ π 2 − K, we find where and E k is the k−th moment of the energy distribution, defined as VIII. APPENDIX: PROOF OF REALIZABILITY OF THE ISOTROPIC ENSEMBLE Theorem: Given any choice of a real vector f 2 orthogonal to f 1 in (34), there exists a set of barriers B ij with relaxation dynamics obeying detailed balance (4) having F −1/2 f 2 as its slowest relaxation eigenvector.
Proof: For our purpose we need to demonstrate at least one choice of barriers. We first note that for any (symmetrized) form of the drivingR with a steady state distribution f 1 , we can obtain, using (4) and (8), formally, a set of barriers as: The only requirement for these B ij to be consistent with our relaxation dynamics is that B ij is a real and symmetric matrix. In other words, it is sufficient thatR ij is symmetric, and thatR ij > 0 for all i = j (note that R ii is then uniquely determined by the condition that f 1 is an eigenvector with eigenvalue 0).
We now show that we can make such a choice for any f 2 . To do so we consider first an initial set of barriers B ij = E i + E j . An explicit calculation shows that the resulting dynamics has a single zero eigenvalue associated with f 1 , and that the rest of the eigenvalues are all −Z(T b ). In this case we haveR ij = e − β b (E i +E j ) 2 . In particular any choice of f 2 orthogonal with f 1 is immediately an eigenvector ofR. It remains to break the degeneracy between f 2 and the other vectors orthogonal to f 1 . We do this by adding a small perturbation toR: This change will only affect the eigenvalue associated with f 2 , changing it to −Z(T b ) + , making it a non degenerate eigenvector.
Clearly, for small enough the positivity ofR ij for i = j will not be affected and the formula (66) will give us a valid set of barriers. (It is enough to take < min ij (e − β b (E i +E j ) 2 )). QED.
Of course, the above procedure yield a very particular type of barriers for each f 2 . There are numerous ways to set up other barriers consistent with a given f 2 .