Out-of-equilibrium clock model at the verge of criticality

TL;DR

This paper studies an out-of-equilibrium 2D clock model that acts as a thermal machine, showing optimal performance near the critical temperature where phase transitions occur, leveraging critical fluctuations to enhance efficiency.

Contribution

It introduces a dynamical response function to quantify the susceptibility of an out-of-equilibrium clock model and demonstrates optimal machine operation at the critical boundary between phases.

Findings

Optimal machine performance occurs near the critical temperature.

Critical fluctuations enhance the system's thermodynamic response.

Out-of-equilibrium driving induces directed rotation of spins.

Abstract

We consider an out-of-equilibrium lattice model consisting of 2D discrete rotators, in contact with heat reservoirs at different temperatures. The equilibrium counterpart of such model, the clock-model, exhibits three phases; a low-temperature ordered phase, a quasi-liquid phase, and a high-temperature disordered phase, with two corresponding phase transitions. In the out-of-equilibrium model the simultaneous breaking of spatial symmetry and thermal equilibrium give rise to directed rotation of the spin variables. In this regime the system behaves as a thermal machine converting heat currents into motion. In order to quantify the susceptibility of the machine to the thermodynamic force driving it out-of-equilibrium, we introduce and study a dynamical response function. We show that the optimal operational regime for such a thermal machine occurs when the out-of-equilibrium disturbance is applied around the critical temperature at the boundary between the first two phases, namely where the system is mostly susceptible to external thermodynamic forces and exhibits a sharper transition. We thus argue that critical fluctuations in a system of interacting motors can be exploited to enhance the machine overall dynamic and thermodynamic performances.

Equations

$$H(\left\{{\theta }_i\right\})=-\frac{k}{2}\sum _{⟨i,j⟩}\cos \left[{\theta }_i-{\theta }_j+\phi \right],$$

$${\chi }_{\omega }={\frac{\partial \omega }{\partial \Delta T}|}_{\Delta T=0},$$

$$⟨m⟩\equiv \frac{⟨M⟩}{L^2}=\frac{⟨|{\sum }_{j=1}^{L^2}\exp (i{\theta }_j)|⟩}{L^2},$$

$$U_L=1-\frac{⟨m^4⟩}{3⟨m^2{⟩}^2}.$$

$$V_L=1-\frac{⟨E^4⟩}{3⟨E^2{⟩}^2},$$

Full text

The critical behaviour of equilibrium systems is well understood and supported by innumerous experimental observations. As a system approaches a critical point larger fluctuations in the local order parameter take longer time to emerge and longer time to decay. Also the response to a change in an external thermodynamic force can become significantly large. For example, in a gas at equilibrium with its liquid phase at the critical point a small change in the pressure produces a large density change, corresponding to a diverging compressibility. Other well known examples of diverging responses are the specific heat and the susceptibility of a critical 3D Ising model. Thus, criticality is characterized by wide fluctuations with diverging characteristic lengths, susceptibilities, and relaxation times.

One may wonder what is the fate of the equilibrium critical fluctuations if a disturbance is applied, that drives the system out of equilibrium. The system must be certainly removed from criticality, but at least for small non-equilibrium perturbation one may expect that fluctuations in the system persist and are still quite large. An even more interesting question is what happens to the diverging response of a critical system when such a non-equilibrium disturbance is applied. Is the response of a system close to its critical points enhanced with respect to the case where the system is non critical?

This question is particularly relevant in the context of stochastic and quantum thermodynamics, where several groups have addressed the issue of performance optimization in microscopic machines operating as engines and motors Verley et al. (2014); Polettini et al. (2015); Polettini and Esposito (2017); Suñé and Imparato (2019). In particular, a few publications have shown that collective effects can enhance the efficiency of interacting microscopic machines. In models of interacting molecular motors the efficiency at maximum power of the many-motor system can be larger than in the single motor case Golubeva and Imparato (2012, 2013, 2014). Models of interacting work-to-work transducers with all-to-all interaction exhibit a similar behaviour Imparato (2015); Herpich et al. (2018); Herpich and Esposito (2019). Furthermore, in ref. Campisi and Fazio (2016) the authors showed that the Carnot efficiency can be approached at finite output power in a system of microscopic interacting Otto engines when the working substances operate at the verge of a second order phase transition and the critical exponents satisfy some given constraints.

In this paper we investigate the role of criticality in the performance of thermal machines driven by temperature gradients: in particular we show that the out-of-equilibrium dynamic response of a system is maximal in the case in which the driving perturbation is applied at the verge of one of its critical points, i.e. when the system is most susceptible to equilibrium thermodynamic forces. We also demonstrate that phase transitions with different strengths are characterized by different dynamical responses, with the smoother equilibrium phase transition exhibiting the feebler out-of-equilibrium response.

Specifically we consider a model of interacting autonomous thermal motors, that converts steady state heat currents into mechanical motion. Autonomous motors are systems which continuously and cyclically convert one form of energy into another: such machines are characterized by the absence of one or more external agents that change their Hamiltonian or temperature. As such they are different from reciprocating engines that perform thermodynamic cycles, such as the Carnot, the Stirling and the Otto cycle, which have been widely investigated in the context of stochastic Seifert (2012) and quantum Kosloff (2013) thermodynamics.

Our working tool will be the 2D clock model, consisting of $L\times L$ discrete rotors on a discrete 2D lattice, interacting through the Hamiltonian

[TABLE]

with $\theta_{i}=2\pi n_{i}/N_{s}$, $n_{i}=0,\dots N_{s}-1$, and with nearest neighbor interactions. We subdivide the lattice into two sublattices, each in contact with a heat bath at a different temperature $T_{\pm}=T_{m}\pm\Delta T$. Specifically we take a 2D square-lattice with periodic boundary conditions and $L$ even, so that each “+” (“-”) rotator has 4 nearest neighbour “-” (“+”) rotators sup . The model eq. (1) is chiral since we take $\varphi>0$ when the spin $i$ belongs to the “-” (cold) sublattice, and $j$ belongs to the ”+” (hot) one.

This model has the advantage that i) a single dimer consisting of only two spins interacting with the Hamiltonian (1) ($i=1,\,j=2)$ has been shown to work as an autonomous motor: when $\varphi\neq n\pi/N_{s}$ and when in contact with two heat baths at different temperatures the spins rotate with a given average frequency, thus converting heat currents into motion and eventually into work Fogedby and Imparato (2017); Hovhannisyan and Imparato (2019); and ii) its non-chiral equilibrium counterpart ($\varphi=0$, $\Delta T=0$) exhibits a rich phase diagram with two second-order phase transitions Lapilli et al. (2006).

In the following we use Monte Carlo (MC) simulations with Metropolis algorithm to first investigate the equilibrium properties of the model (1) ($\Delta T=0$) and then its out-of-equilibrium dynamic properties, in particular the spin rotation rate for $\Delta T\neq 0$, and the heat rates exchanged by the system with the two reservoirs $\dot{Q}_{\pm}$. We define the angular velocity per spin as $\omega=(\omega_{+}+\omega_{-})/2$, where $\omega_{\pm}$ is the angular velocity per spin in the hot (cold) sublattice. The heat currents $\dot{Q}_{\pm}$ along a single MC trajectory can be evaluated by sampling the total energy exchanged by the system with each of the two reservoirs along that trajectory.

We first study the equilibrium behaviour of the model (1), and we benchmark our results with those discussed by Lapilli et al. in Lapilli et al. (2006), where the equilibrium properties of the non–chiral clock model ($\varphi=0$) were extensively studied, and the model phase diagram was established by using MC simulations. In particular, as discussed by the authors, the model (1) with $4<N_{s}<\infty$ and $\varphi=0$ exhibits three phases: a low-temperature ordered and a high-temperature disordered phase, as in the Ising model, along with a quasiliquid intermediate phase. As $N_{s}\to\infty$ the ordered phase is suppressed, and the model becomes equivalent to the continuous $XY$ model. In the following we set the energy scale by taking $k=1$, and we notice that, because of the additional $1/2$ prefactor in our model (1), the temperature scale in the present paper is halved with respect to Lapilli et al. (2006). We label $T_{1}^{\mathrm{eq}}$ the equilibrium transition temperature between the low-temperature ordered and the quasi-liquid phase, and $T_{2}^{\mathrm{eq}}$ the equilibrium transition temperature between the quasi-liquid and the disordered phase. By using MC simulations, Lapilli et al. Lapilli et al. (2006) provided an expression for the transition temperature $T_{1}^{\mathrm{eq}}$ as a function of the number of states for $N_{s}\geq 6$. They also reported an extended universality, which implies that the transition $T_{2}^{\mathrm{eq}}$ must be of the Berezinskii–Kosterlitz–Thouless (BKT) class, thus indistinguishable from $N_{s}=\infty$, for $N_{s}\geq 8$. Such equilibrium temperatures are reported in fig. 1 as lines, together with the results from our MC simulations.

Similarly, we run equilibrium MC simulations ($\Delta T=0$) for the chiral model $\varphi=\bar{\varphi}$, where we fix $\bar{\varphi}=\pi/(2N_{s})$, and evaluate the equilibrium transition temperatures $T_{1}^{\mathrm{eq}}$ and $T_{2}^{\mathrm{eq}}$, for different values of $N_{s}$. This is done through finite-size scaling analysis of the order parameter $\langle m\rangle\equiv\langle|\sum_{l}\exp({i\theta_{l}})|\rangle/{L^{2}}$ (the sum runs over all the rotators in the system), see sup ; Lapilli et al. (2006) for further details. The results for the two temperatures are also shown in fig. 1 and the agreement with the findings of Lapilli et al. (2006) are excellent for $N_{s}\geq 8$. This can be understood by noticing that the chosen phase $\varphi=\bar{\varphi}$ decreases with $N_{s}$, so the equilibrium properties of the two models become equivalent in the limit of large $N_{s}$.

By using finite-size scaling analysis, we find that the critical exponent for the order parameter at the transition temperature $T_{1}^{\mathrm{eq}}$ is a decreasing function of the number of states, thus signaling a sharper transition at $T_{1}^{\mathrm{eq}}$ as $N_{s}$ increases, see Supplemental material sup .

We now turn our attention to the model out-of-equilibrium properties. As for the case of two isolated spins interacting with the Hamiltonian (1) discussed in Hovhannisyan and Imparato (2019), also for the present model with $L^{2}$ spins there is a broken rotational symmetry when $\varphi\neq n\pi/N_{s}$, see sup .

Thus, as in the 2 spin model Fogedby and Imparato (2017); Hovhannisyan and Imparato (2019), when the broken rotational symmetry is combined with a constant temperature gradient, the system is driven into an out-of-equilibrium steady state, characterized by non-vanishing rotation rate per spin $\omega$. The average angular velocity of the out-of-equilibrium system with $T_{\pm}=T_{m}\pm\Delta T$ is plotted in fig. 2 as a function of $T_{m}$ for $N_{s}=16$ and for different values of $L$ and $\Delta T$. We find that at large $L$ the curve $\omega(T_{m})$ becomes independent of $L$, see fig. 2-(a).

In order to investigate the role of the critical fluctuations in a non-equilibrium system we introduce a dynamic response to a non-equilibrium perturbation, akin to the susceptibility in equilibrium systems. In equilibrium statistical mechanics the susceptibility of an order parameter $\phi$ with respect to its conjugate thermodynamic force $f_{\phi}$ reads $\chi_{\phi}=\partial\phi/\partial f_{\phi}|_{f_{\phi}=0}$. Similarly we introduce the susceptibility

[TABLE]

i.e., the response of the angular velocity with respect to the thermodynamic force that generates it. In the following we argue that this is the relevant quantity to characterize the out-of-equilibrium response of a system close to its critical points.

Given that $\omega=0$ when $\Delta T=0$ we can thus evaluate the response function (2) by numerical derivation $\chi_{\omega}\simeq\omega/\Delta T$. Unless otherwise stated, in the following we take $\Delta T\simeq T_{1}^{\mathrm{eq}}/100$, for any $N_{s}$. Indeed when $\Delta T$ is as small as $T_{1}^{\mathrm{eq}}/100$ the curves $\omega(T_{m})/\Delta T$ collapse into a single curve, see fig. 2-(b).

The results for $\chi_{\omega}$ as a function of $T_{m}$ and for different values of $N_{s}$ are shown in fig. 3-(a). We find that each curve $\chi_{\omega}(T_{m},N_{s})$ exhibits a peak at a temperature that depends on $N_{s}$. We indicate such a temperature $T^{*}$: it decreases as $N_{s}$ increases. We refer to the value of the dynamic response function at its maximum as $\chi_{\omega}^{*}=\chi_{\omega}(T^{*},N_{s})$.

For any number of states $N_{s}$, we then compare the temperature $T^{*}$, as obtained for two values of $\Delta T$, with the model equilibrium properties. In particular we compare the peak temperature with the two transition temperatures $T_{1}^{\mathrm{eq}}$ and $T_{2}^{\mathrm{eq}}$. Interestingly, the temperature $T^{*}$ (triangles in fig. 1) where $\chi_{\omega}(T_{m},N_{s})$ has its maximum $\chi_{\omega}^{*}$ follows very closely the equilibrium transition temperature $T_{1}^{\mathrm{eq}}$ for $N_{s}\geq 4$ (full circles in fig. 1) with a relative difference $(T^{*}-T_{1}^{\mathrm{eq}})/T_{1}^{\mathrm{eq}}\lesssim 0.15$. This in turn implies that the points $T^{*}$ with $N_{s}\geq 6$ lie near the curve given in Lapilli et al. (2006) for the equilibrium temperature $T_{1}^{\mathrm{eq}}$ of the non-chiral model (full line in fig. 1). Accordingly, we conclude that the response of the system to the out-of-equilibrium perturbation $\Delta T$ is maximal when the perturbation is applied on the system in the proximity of one of its critical points: namely the transition point between the ordered and the quasiliquid phase. A more precise analysis on the correspondence between $T^{*}$ and $T_{1}^{\mathrm{eq}}$ can be found in sup .

Unlike $T_{1}^{\mathrm{eq}}$, the temperature $T_{2}^{\mathrm{eq}}$ does not give rise to any nonvanishing response, for any value of $N_{s}$. This might be seen as a consequence of the fact that the critical behaviour of the $XY$ model at the BKT transition exhibits very weak singularities Kosterlitz and Thouless (1973); Chaikin and Lubensky (1995). For example the specific heat at $T_{\mathrm{BKT}}$ displays a singularity which is essentially unobservable Chaikin and Lubensky (1995).

In order to quantify the performance of the system, and its capability of transforming heat into rotational motion, we introduce the figure of merit $z=\langle\omega\rangle/\langle\dot{Q}_{+}\rangle$, where $\langle\dot{Q}_{+}\rangle$ is the steady state heat current entering the system from the hot reservoir. Interestingly, such figure of merit follows the same behaviour of $\omega$ and $\chi_{\omega}$, that is, it peaks around $T^{*}$. This result confirms that the optimal operational regime is achieved at $T_{m}\simeq T_{1}^{\mathrm{eq}}$, see fig. 3-(b). This behaviour can be understood by noticing that while the angular velocity exhibits a somewhat abrupt peak around $T^{*}$ (figs. 2 and 3-(a)), the heat current $\langle\dot{Q}_{+}\rangle$ presents a smoother increase and decrease in that temperature range (data shown in sup ). Thus their ratio $z$ exhibits the same peaking behaviour as $\omega$ and $\chi_{\omega}$.

Inspection of fig. 3 suggests that the height of the peak in both the dynamic susceptibility and the figure of merit is an increasing function of the number of states $N_{s}$. Thus as next point we investigate the dependency of $\chi_{\omega}^{*}$ on $N_{s}$; the results are plotted in fig. 4. Such a quantity increases linearly with the number of states, although with a slow rate. However, it is worth noticing that the velocity goes to zero as $N_{s}\to\infty$, as the rotational symmetry is restored (see sup ), as can be seen in the lower inset of fig. 4. Therefore the continuous model $N_{s}\to\infty$ exhibits a diverging $\chi^{*}_{\omega}$, while $T^{*}\to 0$.

Interestingly, the increase of $\chi_{\omega}^{*}$ is accompanied by the decrease of the order parameter critical exponent for the corresponding equilibrium system at $T_{1}^{\mathrm{eq}}$, see  sup . We thus conclude that as the equilibrium transition becomes sharper, so does the corresponding out-of-equilibrium response.

To summarize, we have used the 2D chiral clock model to show that the optimal operational regime of a system of interacting autonomous motors occurs at the verge of a phase transition, in the limit of small driving thermodynamic forces. While the out-of-equilibrium response around the ordered-quasiliquid phase transition is noticeably sharp, the system response is vanishing in the other temperature ranges, in particular at the quasiliquid-disordered phase transition temperature. The same behaviour is found if one considers a figure of merit akin to a miles per gallon equivalent: the motors’ “consumption” is most favorable around the strong phase transition. Interestingly, we find also that the response of the system increases, although slowly, with the number of states (directions) $N_{s}$ in the clock model. This increase occurs while the corresponding phase transition for the equilibrium system becomes sharper. Our results provide strong evidence for the protocol that one has to follow in order to find the optimal working regime in systems of microscopic interacting motors: one has to first identify the system critical points, if any, and then systematically perturb the system at the verge of criticality.

Experimentally, a 2D system of interacting molecular rotors deposited on an hexagonal lattice have been reported in Zhang et al. (2016). That artificial molecular system represents a possible setup for an experimental realization of the model we have studied here.

Appendix A Supplemental Material: Out-of-equilibrium clock model at the verge of criticality

A.1 Critical points for the chiral clock model at equilibrium

The Hamiltonian for the clock model in two dimensions eq. (1) with $\varphi=0$ has a global symmetry $Z_{N_{s}}$ ($H(\{\theta_{i}\})$ is invariant under the transformations $\theta_{i}\to\theta_{i}+2\pi n/N_{s}$). Thus one introduces the order parameter

[TABLE]

that is invariant under the same symmetry group as the Hamiltonian in the high-temperature disordered phase. Besides this phase, the system exhibits two other phases for $N_{s}\geq 4$ (low-temperature ordered and quasi-liquid intermediate phase) and thus two critical points, as discussed in Lapilli et al. (2006): an ordered to quasi-liquid second-order phase transition and a quasi-liquid to disordered phase transition, that turns out to be of the Berezinskii-Kosterlitz-Thouless class for $N_{s}\geq 8$.

In the present paper we have considered the chiral model with $\varphi=\bar{\varphi}$. Similarly to the model with only two rotators discussed Hovhannisyan and Imparato (2018), also for the present model with $L^{2}$ rotors there is a broken rotational symmetry, in the sense that when $\varphi\neq n\pi/N_{s}$, one cannot find a discrete sublattice rotation $l^{+}$ (or $l^{-}$) such that $H(-\bm{\theta})=H(\{\theta^{+}_{i}+2\pi l^{+}/N_{s},\theta^{-}_{j}\})$, $\forall\bm{\theta}$ (or $H(-\bm{\theta})=H(\{\theta^{+}_{i},\theta^{-}_{j}+2\pi l^{-}/N_{s}\})$), where $\bm{\theta}=\{\bm{\theta}^{+},\,\bm{\theta}^{-}\}$ is an arbitrary configuration of the system, and $\bm{\theta}^{+},\,\bm{\theta}^{-}$ are the configurations of the two sublattices, see fig. 5. This symmetry plays a fundamental role in the out-of-equilibrium model described in the main text. In the limit of $N_{s}\to\infty$ the model becomes the continuous $XY$ model and the energy eq. (1) exhibits an $O_{2}$ symmetry, the group of rotations in a two-dimensional plane, and thus the model does not break the spatial symmetry discussed above. The consequence is that the angular velocity for the out-of-equlibrium model vanishes as $N_{s}\to\infty$ as discussed in the main text.

Through equilibrium MC simulations of the model we have evaluated the order parameter (3) as a function of the temperature. Such a quantity is shown in fig. 6 for $N_{s}=6$. Inspection of this figure indicates that the chiral model also exhibits two phase transitions. The order parameter for the standard clock model ($\varphi=0$) is included for comparison.

The order parameter curves display two abrupt drops as the temperature is increased, from the low-temperature ordered phase ($\langle m\rangle\approx 1$) to the high-temperature disordered phase ($\langle m\rangle\approx 0$). In order to evaluate the corresponding transition temperatures we follow the approach used in  Lapilli et al. (2006). To obtain the low-temperature transition $T_{1}^{eq}$, we evaluate the minimum of the temperature derivative of $m$ for different system sizes $L$, $T_{1}^{eq}(L)$, see inset in fig. 6. We then use finite-size scaling analysis Plischke and Bergersen (2006) to obtain $T_{1}^{eq}=T_{1}^{eq}(L\to\infty)$. For the high-temperature transition $T_{2}^{eq}$ we exploit the Binder’s fourth order cumulant Binder and Heermann (2010) of the order parameter,

[TABLE]

For $T>T_{2}^{eq}$ and assuming $L\gg\xi$, where $\xi$ is the correlation length of the order parameter fluctuations, $U_{L}\propto L^{-d}$ ($d$ is the dimension of the system). On the other hand, $U_{L}\to U_{\infty}=2/3$ when $T<T_{2}^{eq}$ and $L\gg\xi$; and $U_{L}$ varies weakly with $T$ and $L$ for $L\ll\xi$. Thus, the cumulant (4) is scale invariant at the transition temperature, and therefore the transition temperature is given by the intersection between ratios of cumulants for different system sizes Binder and Heermann (2010).

The critical temperatures as computed employing the aforementioned approaches for the equilibrium chiral clock model are depicted in fig. 1. The second equilibrium transition temperature is the same for both the chiral version ($\varphi\neq 0$) and the standard clock model ($\varphi=0$), as it can be anticipated by the collapse of the corresponding order parameters curves in fig. 6. However, the order parameter curves deviate as they approach the lower-temperature critical point, thus the chiral clock model exhibits the first transition at a lower temperature. This deviation decreases with $N_{s}$ as the two models become equivalent in this limit.

For different values $N_{s}$, and different sizes $L$, we evaluate the Binder’s fourth order cumulant for the energy Binder and Heermann (2010)

[TABLE]

around the first transition temperature $T_{1}^{\mathrm{eq}}$. The latent heat is proportional to $\lim_{L\to\infty}(2/3-\min_{T}V_{L})$ as discussed in Binder and Heermann (2010). We find that $\lim_{L\to\infty}\min_{T}V_{L}=2/3$, which thus signals a second-order phase transition. We then proceed to evaluate the critical exponent $\beta$ for the order parameter $(m-m_{c})\propto(T-T_{1}^{\mathrm{eq}})^{\beta}$, where $m_{c}$ is the value of $m$ at the transition point. We use finite-size scaling analysis of the order parameter for different sizes $L$ Plischke and Bergersen (2006) to obtain the critical exponent $\beta$. The results are shown in fig. 7, where we plot the exponent $\beta$ as a function of $N_{s}$. We find that $\beta$ is a decreasing function of the number of states $N_{s}$, indicating that the transition becomes sharper as the number of states increases. This result has to be compared with the results for the dynamical susceptibility in the main text, whose maximum increases with $N_{s}$, see in particular fig. 4.

A.2 Correspondence between $T^{*}$ and $T_{1}^{\mathrm{eq}}$

The temperature at which the susceptibility eq. (2) is maximal $T^{*}$ has been shown to depend on $\Delta T$ , see fig. 2-(b). Its behavior when varying the number of states $N_{s}$ is depicted in fig. 1, where two different temperature gradients are considered $\Delta T\simeq T_{1}^{\mathrm{eq}}/10$ and $\Delta T\simeq T_{1}^{\mathrm{eq}}/100$. The maximum temperature $T^{*}$ follows closely the equilibrium clock model transition temperature $T_{1}^{\mathrm{eq}}$ for both $\varphi=0,\,\bar{\varphi}$), yet some discrepancy is already anticipated in fig. 1. In order to compare $T^{*}$ and $T_{1}^{\mathrm{eq}}$ a more convenient representation of the data in fig. 1 is to plot the relative deviation of the maximum susceptibility’s temperature $T^{*}$ with respect to the equilibrium transition temperature $T_{1}^{\mathrm{eq}}$ as a function of the number of states, see fig. 8.

When decreasing the number of states, and hence when the standard clock model begins to deviate from the chiral one, the maximal susceptibility temperature $T^{*}$ tends to approach the corresponding equilibrium temperature $T_{1}^{\mathrm{eq}}$. Yet a deviation persists, so that $T^{*}$ is always above $T_{1}^{\mathrm{eq}}$ when the temperature gradient is small. On the other hand, the difference between $T^{*}$ and $T_{1}^{\mathrm{eq}}$ is larger for the smaller $\Delta T$ considered, in agreement with the results plotted in fig. 2-(b).

A.3 Heat current

The heat currents $\dot{Q}_{\pm}$ along a single MC trajectory can be retrieved by exploiting the conservation of energy: the energy rate exchanged by each of the two sublattices with the corresponding heat reservoir must amount to the energy variation due to the sublattice’s spin transitions, given by eq. (1). Accordingly, the sum of the two heat currents should sum up to zero on average.

The heat current as a function of $T_{m}$ for a given value of $N_{s}$ reaches a maximum which depends on $N_{s}$, similarly to what happens for the angular velocity, see fig. 9-(a). Yet a close inspection shows that the curve reaches such a maximum at a larger temperature than $T^{*}$, see fig. 9-(b), where we compare the angular velocity $\langle\omega\rangle$ and the heat current $\langle\dot{Q}_{+}\rangle$ for two values of $N_{s}$. The angular velocity exhibits a sudden rise/decrease around $T^{*}$. The heat current has a somewhat smoother rising and decreasing behaviour around the maximum. We see that the system still transports heat ($\langle\dot{Q}_{+}\rangle\neq 0$) in the range of $T_{m}$ where $\langle\omega\rangle$ vanishes due to the wandering of the rotators between their different positions, although this movement does not result in a net rotation. The combination of the abrupt variation of $\langle\omega\rangle$ around $T^{*}$ together with the smoother behaviour displayed by $\langle\dot{Q}_{+}\rangle$ implies that the figure of merit $z$ (i.e. their ratio) follows the same behaviour of the angular velocity and peaks around $T^{*}$, see fig. 9-(c).

The heat current’s peak at $T\gtrsim T^{*}$ is nonetheless not unique; a second maximum (possibly larger than the one discussed above) is found above $T_{2}^{\mathrm{eq}}$, see fig. 10-(a), although only for small $N_{s}$. The maxima displayed by the curves $\langle\dot{Q}_{+}\rangle/\Delta T$ close to the transition temperature $T_{1}^{\mathrm{eq}}$ is constant with respect to the number of states $N_{s}$, see fig. 10-(a), at variance with the mechanical current response function as defined in eq. 2, see also fig. 3-(a) in the main text. A rescaling of the x-axis thus leads to a collapse of the heat current curves around the transition temperature $T_{1}^{\mathrm{eq}}$, see fig. 10-(b). Both the mechanical and thermal currents are expected to be linear with respect to the thermodynamic force when $\Delta T$ is small, as it is in our case, according to the Onsager reciprocal relations. However, we find a different scaling behavior of the two currents with respect to $N_{s}$, as suggested when comparing fig. 3-(a) and fig. 10-(a).