Entropy production at criticality in a nonequilibrium Potts model

Understanding nonequilibrium systems and the consequences of irreversibility for the system's behavior as compared to the equilibrium case, is a fundamental question in statistical physics. Here, we investigate two types of nonequilbrium phase transitions, a second-order and an infinite-order phase transition, in a prototypical q-state vector Potts model which is driven out of equilibrium by coupling the spins to heat baths at two different temperatures. We discuss the behavior of the quantities that are typically considered in the vicinity of (equilibrium) phase transitions, like the specific heat, and moreover investigate the behavior of the entropy production (EP), which directly quantifies the irreversibility of the process. For the second-order phase transition, we show that the universality class remains the same as in equilibrium. Further, the derivative of the EP rate with respect to the temperature diverges with a power-law at the critical point, but displays a non-universal critical exponent, which depends on the temperature difference, i.e., the strength of the driving. For the infinite-order transition, the derivative of the EP exhibits a maximum in the disordered phase, similar to the specific heat. However, in contrast to the specific heat, whose maximum is independent of the strength of the driving, the maximum of the derivative of the EP grows with increasing temperature difference. We also consider entropy fluctuations and find that their skewness increases with the driving strength, in both cases, in the vicinity of the second-order transition, as well as around the infinite-order transition.


Introduction
Phase transitions are ubiquitous in nature and generally occur in equilibrium as well as nonequilibrium systems. In either case, the transition is due to internal interactions and often goes along with the breaking of spatial symmetries as a reaction to the variation of a control parameter below its critical value, detectable by the emergence of an appropriate order parameter. Phenomena like the occurrence of phase transitions in one spatial dimension are solely observed in nonequilibrium systems [1][2][3][4][5][6][7][8][9][10][11][12][13]. In contrast, other properties related to phase transitions are identical and thus do not allow perceiving whether a system is in a state of thermal equilibrium or not. In equilibrium, it is well established that continuous (second-order) phase transitions are accompanied by power-law divergences of multiple measurable quantities, such as the magnetic susceptibility or the spin-spin correlation length (for Ising-like models) and there are already numerous examples for nonequilibrium systems that can be characterized in this manner as well [14][15][16][17][18][19][20]. However, a general theory for nonequilibrium phase transitions is still missing, and it is not per se clear whether the critical exponents of a system stay the same (i.e., the system remains in the same universality class) when driven away from equilibrium.
To analyze nonequilibrium systems there exists another, yet almost completely unconnected, tool, that is, the entropy production (EP). This quantity is strictly positive for nonequilibrium and exactly zero for equilibrium systems. The total EP is a fundamental quantity of great importance in statistical physics that already plays a central role in (stochastic) thermodynamics and information theory, and it is known to fulfill various laws, including the famous fluctuation theorems [21][22][23][24][25][26][27]. Moreover, the EP can be defined in a very general manner (not system-specific) and is, in principle, a meaningful quantity for any complex system. This is because it solely depends on (state and transition) probabilities, and does not rely on concepts like energy or temperature. Thus, it can be defined and calculated also, e.g., in non-physical systems, like social dynamics, or opinion formation, which are at the same time known to undergo phase transitions [28,29]. For these reasons it is tempting to investigate this quantity with regard to critical behavior in nonequilibrium systems.
A few recent studies have already started to investigate the behavior of total EP around criticality in different lattice-based models by means of Monte-Carlo simulations and mean field theory [30,31]. For example, for a one-dimensional KPZ interface growth model [32,33] the rate of EP around the critical point of a first-order phase transition was calculated. Furthermore, there are several studies considering spin systems with up-down (Z 2 ) symmetry. These include an interacting lattice gas model in contact with two heat and particle reservoirs [30], the majority vote model [34,35], and an Ising model externally driven by an deterministically oscillating magnetic field [36]. In all cases, the EP rate was either found to jump or display an inflection point at criticality, such that its first derivative with respect to an appropriate control parameter exhibits marked behavior around the critical point. In particular, the derivative shows a discontinuity or a power-law divergence, bearing a resemblance with the susceptibility, the spin-spin correlation length and the specific heat. Of special interest for the present work is a study considering a variant of the square lattice Ising model with nearest-neighbor interactions, whose spins are in contact with two heat baths at temperatures T 1 and T 2 in a checkerboard arrangement [37,38]. In this study, the derivative of the EP was found to diverge around the second-order phase transition with the same critical exponent as the specific heat [38]. This raises the question whether these diverging quantities generally behave alike at criticality.
In the present paper, we aim to generalize the previous findings to spin symmetries different from Z 2 , and to other types of phase transitions. To this end, we consider a nonequilibrium q-state vector Potts model around criticality. Depending on the value of q, the model displays either a second-order phase transition from a paramagnetic (PM) to a ferromagnetic (FM) phase or an infinite-order phase transition [39][40][41] from a PM to a quasi long-range ordered Berezinskii-Kosterlitz-Thouless (BKT) phase. We investigate the model in the vicinity of both types of phase transitions under nonequilibrium conditions. In particular, we couple the spins to heat baths at two different temperatures T 1 and T 2 in such a way, that all nearest-neighbors of each spin are in contact only with spins coupled to the respective other heat bath. Hence, two sublattices are formed in a checkerboard arrangement similar to [37,38]. As a consequence of the two involved temperatures, a net heat flow from the hotter to the colder heat reservoir is induced. This setup drives the system in a nonequilibrium steady state where it constantly produces entropy, which is exported to the environment. We investigate the system numerically using Monte Carlo simulations with Glauber dynamics. To deepen our understanding of nonequilibrium phase transtions, we study the EP and its fluctuations, as well as standard quantities such as the magnetization. Using the finite-size scaling technique [42][43][44][45] to carefully analyze the critical behavior based on numerical data, we dedicate a detailed analysis to the question whether the specific heat and the derivative of the EP behave alike at criticality for both types of phase transitions.
Our analysis reveals that the derivative of the EP rate with respect to temperature in the PM disordered phase of the 4-state vector Potts model shares indeed similarities with the specific heat. Specifically, it shows power-law behavior as function of the distance from criticality, however, in contrast to the specific heat, with non-universal scaling exponent. Additionally, its maximum value diverges at criticality and also shows power-law behavior as function of system size with a non-universal scaling exponent. In contrast, for the XY model with q → ∞, the EP rate resembles the behavior of the specific heat. Both quantities do not diverge in the vicinity of the transition from the PM to the BKT phase.

The q-state vector Potts model
The Hamiltonian of the q-state vector Potts model (or the q-state clock model) with nearest-neighbor spin interactions on a discrete lattice without any externally applied magnetic field is defined as Here, J represents the coupling constant between interacting spins which is set to unity and thereby favours ferromagnetic order. The sum in Eq. (1) runs over all neighboring lattice sites ij . We consider square lattices in two dimensions with a total number of L 2 spins, where L the lateral extension of the system which we refer to as the "system size". The spins σ i = e x cos(θ i ) + e y sin(θ i ) = [cos(θ i ), sin(θ i )] are represented as two-component unit vectors in the x − y plane located on discrete Figure 1. Illustration of the q-state vector Potts model on a two-dimensional square lattice with periodic boundary conditions. The spins are coupled to heat baths at two different temperatures T 1 and T 2 , here indicated by the colors red and blue, respectively, with a checkerboard arrangement. Black lines represent the nearest-neighbor interactions of strength J, while red and blue lines highlight the two sublattices formed by the coupling to heat baths of different temperatures. and equidistant positions i on the lattice. Within the q-state vector Potts model, the possible angles θ i ∈ [0, 2π] of the spins are given by where the integers a = 0, 1, 2, ..., q − 1 determine the possible orientations of the spins. At q = 2, the model reduces to the classical Ising model with up-down (Z 2 ) spin symmetry, whereas in the limiting case q → ∞, the model corresponds to the XY model where the spin orientations are continuous within the plane. In what follows, we focus on the cases q = 4 and q → ∞. For, q = 4 the model (i.e., the Ashkin-Teller model) shows a second-order phase transition from a PM to a FM phase similar to the Ising model, yet with different characteristics, i.e., different critical exponents [46]. In contrast, in the two-dimensional XY model (where q → ∞), there exists no long-range ordered FM phase at finite temperatures as stated by the Mermin-Wagner theorem [47]. Instead, the sytem undergoes an infinite-order BKT transition from a PM to a BKT phase.
In order to study the dynamical evolution of the system in presence of thermal noise, we perform Monte-Carlo simulations with single spin-flip Glauber dynamics. Here, the rate w µν (i) for a transition of a randomly chosen spin σ i from state ν (before the spin flip) to µ (after the spin flip) depends only on the temperature T k of the heat bath the considered spin is coupled to, as well as on the energy difference ∆H = H(ν) − H(µ) related to a flip of spin σ i , that is, If all spins σ i are exposed to a single heat bath at temperature T , the system (that is initially prepared in a configuration with random spin orientations) eventually reaches a state of thermal equilibrium and thus, does not produce entropy, Π = 0 (see Sec. 3 for a definition of entropy production). In contrast, here we drive the system into a nonequilibrium steady state by coupling the spins σ i to two different heat baths T k (k = 1, 2) which are kept at temperatures T 1 and T 2 . With this setup, the system is out of equilibrium whenever T 1 = T 2 . There is a constant heat fluxQ from the hotter to the colder heat reservoir that goes along with a constant rate of entropy production, Π > 0 (see Sec. 3). Our setup splits the system into two sublattices, L 1 and L 2 , each containing all spins connected to the bath at T 1 or T 2 , respectively. All four nearest-neighbors of a spin σ i are coupled to the respective other heat bath, yielding a checkerboard configuration as illustrated in Fig. 1. In the following, we fix T 2 , but vary T 1 and calculate all observables as function of the mean temperature T = (T 1 + T 2 )/2.

Entropy production in the vector Potts model
A key quantity that distinguishes systems out of thermal equilibrium from those in equilibrium is the constant net production of entropy. In general, the dynamical evolution of physical systems that possess a finite set ν ∈ Ω of discrete microstates (e.g., the spin configurations of the vector Potts model at finite q) can be described as continuous time Markov Chains. For such systems, the time-dependent system (or Shannon) entropy [48] is given by the Boltzmann-Gibbs expression Here, the sum runs over all microstates (i.e., spin configurations) and p ν (t) represents the occupation probability of state ν at time t. By coupling the system to an infinitely large heat bath at temperature T (that always maintains a state of thermal equilibrium), we can formulate an expression for the time-dependent change of entropy [24] that one the one hand, originates from the total system internal production of entropy Π(t) and, on the other hand, by the exchange of entropy Φ(t) with the environment, In order to formulate explicit expressions for Π(t) and Φ(t), we make use of the fact that the (generally time-dependent) occupation probabilities p ν (t) obey a master equation The change ∂ t p ν (t) stems first, from the total incoming probability flow µ =ν w νµ (t)p µ (t) consisting of all possible state transitions µ → ν happening with transition rates w νµ (t). The second contributive to ∂ t p ν (t) is the total outgoing probability flow µ w µν (t)p ν (t) due to transitions ν → µ.
For systems in thermal equilibrium, the detailed balance (DB) condition, w νµ p µ = w µν p ν , holds for all µ and ν. When DB is violated, there are non-vanishing local probability flows between certain microstates, i.e., w νµ p µ − w µν p ν = 0. As a consequence, the system constantly produces entropy Π(t) > 0, which is given by [49] Equation (7) obeys the thermodynamically expected properties: Π(t) nullifies in thermal equilibrium, and is strictly positive otherwise, in accordance with the second law of thermodynamics. For nonequilibrium stationary states (no time-dependency), Π(t) = Π, one has ∂ t S = 0 and consequently from Eq. (5), Π = Φ. Additionally, Eq. (7) reduces to Equation (8) can be computed numerically by averaging over many transitions µ → ν from the current state µ of the Markov Chain (i.e., the current spin configuration of the lattice) in the steady state. In spin systems with discrete spin orientations like the vector Potts model with finite q, state transitions µ → ν correspond to the flipping of a randomly chosen spin σ i on lattice site i. Thus, the sum in Eq. (8) can be written as an average over all lattice sites Here, w i νµ [see Eq. (3)] corresponds to the Glauber-type flipping rate of the spin σ i on lattice site i (which is connected to the heat bath at T k ), inducing a transition from the current state µ to state ν due to a change of the current orientation θ i of spin σ i to any other allowed one. The steady exchange of entropy with the environment results from the net heat fluxQ from the hotter to the colder sublattice. We here employ the sign conventionQ > 0 for the heat flow from hot to cold. Due to energy conservation (and because no external fields, forces or further gradients act on the system), all of the three relevant heat flows transport the same amount of energy per timestep: the flow from the hotter heat bath T 1 (here we assume for a moment T 1 > T 2 ) to the corresponding sublattice L 1 , the flow from L 1 to L 2 , and the heat flow from L 2 to the cold bath at T 2 (or everything reversed, if T 2 > T 1 ). This amounts to an overall entropy flow to the environment of Equation (9) can be used to calculate the mean of the entropy production rate, Π, for systems with a finite set Ω of discrete microstates (i.e., the vector Potts model with q = 4) by means of Monte-Carlo simulations. In Sec. 5, we show results for the entropy production rate per spin, which is given by π = Π/L 2 .
A problem with Eq. (9) is that it can not be used to calculate Π in systems with continuous degrees of freedom, as it is the case for the vector Potts model with q → ∞ (XY model). As an alternative, we calculate the entropy production rate by following individual stochastic trajectories consisting of consecutively executed state transitions s n−1 → s n between microstates as the system dynamically evolves via the reorientation of single spins σ i [24]. Such trajectories correspond to a sequence (s 0 → s 1 → s 2 ... → s n−1 → s n ... → s l−1 → s l ) consisting of l state transitions, each connected with a transition rate, w i snsn−1 , given by Eq. (3). Note that each individual state that is part of the trajectory simply corresponds to one of the microstates of the system, s n ∈ Ω, i.e., to one of the possible spin configurations. Since there are infinitely many states for q → ∞, we make use of the fact that each transition s n−1 → s n [for which we know the transition rate w i snsn−1 according to Eq. (3)] along the stochastic path due to a random reorientation of a spin σ i is associated with a (stochastic) entropy exchange with the environment ∆φ = ln w i snsn−1 /w i sn−1sn [24] and an associated local heat exchange of ∆φ/T k . As a consequence, the total change of entropy along an individual stochastic path consisting of l transitions reads In the limit of infinitely long trajectories, l → ∞, Eq. (10) divided by the length l of the trajectory becomes identical to the ensemble averaged medium entropy production rate Φ due to the ergodicity of the system. In a steady state this is further identical to the total entropy production rate as calculated according to Eq. (9) In addition, we also use Eq. (10) to obtain distributions P [∆φ(l = 100)] for the 4-state model as well as the version with q → ∞ (see Sec. 4 for details).

Measurement details and parameter settings
In the present study, simulations of the vector Potts model with nearest-neighbor interactions are performed on square lattices with lateral extension ranging from L = 16 to L = 96. This means that we consider L 2 = 256 up to L 2 = 9216 spins. Before calculating any physical quantity, we first let the system evolve for 5 × 10 4 Monte Carlo steps (MCS) [where one MCS consists of L 2 spin flip attempts with spin-flip rates w i µν according to Eq. (3)] to assure that the system has reached a steady state. We then let each system further evolve up to a maximum of 10 6 MCS and use (depending on the system size L) between 100 and 1000 realizations for each parameter setting (i.e., combination of T 1 and T 2 ) in order to guarantee the convergence of average quantities.
In order to quantify the phase behavior of the system, we calculate the magnetic order parameter The value of m ∈ [0, 1] is a measure for the spin ordering in the system. For perfect order, m = 1, while in a completely disordered system, m = 0. We also define the magnetic order parameters for the two sublattices L 1 and L 2 where the index k = 1, 2 denotes the respective sublattice.
To precisely determine the value of the critical temperature T c where the phase transition (from the PM to FM or to the BKT phase) sets in, we compute the fourthorder Binder cumulant [42] of the magnetic order parameter m which is universal at criticality. The critical value T c of the control parameter is given by the intersection point of U 4 for different lateral sizes L of the system.
In addition, we calculate the specific heat per lattice site which is given by what can also be expresses as d E /dT . Here, E corresponds to the average energy per spin. The specific heat is known to exhibit power-law scaling as the critical critical temperature T c is approached from the PM disordered phase. Additionally, C v peaks at T c where its maximum shows power-law scaling as function of L which is often universal [50]. The average EP rate per spin, π = Π/L 2 , is calculated according to Eq. (10) for the 4-state vector Potts model, while Eq. (11) is used in the case of the XY model (q → ∞) where the spin orientation is continuous. However, π can also be obtained from Eq. (11) for the 4-state model. In both cases, we check whether the change of π, with respect to the control parameter T , dπ/dT , shows universal features (similar to the specific heat C v ) regarding its scaling behavior as function of system size L around the critical point T c of the phase transition.
Distributions P (φ) of the change of entropy φ = ∆φ(l) for trajectories of length l = 100 are obtained via Eq. (10) for q = 4 and q → ∞. We calculate φ for the whole lattice and the individual sublattices L 1 and L 2 , respectively. This is done by defining trajectories that only account for the change of entropy induced by state transitions due to a reorientation of spins which are connected to the respective heat bath T k (k = 1,2). The distributions P (φ) are obtained from at least 10 7 trajectories.

Results
In this section, we present a numerical investigation of the nonequilibrium vector Potts model with discrete (q = 4) and continuous (q → ∞) symmetry. In both cases, we find that the nonequilibrium model exhibits the same type of phase transition as in the equilibrium case. We therefore focus for q = 4 on the transition from the spindisordered PM to the spin-ordered FM phase, whereas for q → ∞ we analyze the BKTlike transition from the disordered PM to the quasi long-range ordered BKT phase. In both cases the transitions are continuous in the order parameter. To characterize the critical behavior, we study the specific heat and the entropy production, and we compare the results in the vicinity of the critical point for both kinds of phase transition.

Phase transition in the discrete vector Potts model with q = 4
Before we numerically investigate the phase transition of our nonequilibrium model, let us briefly review some important properties of the equilibrium version of the 4-state vector Potts model. It is well known [51][52][53] that the equilibrium model exhibits a second-order phase transition at T eq c = 1.13, which is half the critical temperature of the classical Ising model (T eq c = 2.26) that exhibits up-down symmetry. The critical exponents of the 4-state version are different from the Ising model [46].
To begin with the analysis of the nonequilibrium model, we consider the behavior of the ensemble-averaged magnetization m [see Eq. (12)], which serves as a global order parameter. Figure 2 displays m as function of the mean temperature T for four different (fixed) values of T 2 and various system sizes ranging from L = 16 to L = 96. The most prominent observation is that while decreasing the mean temperature from high values, the order parameter increases and eventually approaches its maximum value, m = 1 (reflecting perfect spin order). This implies the existence of a stable FM phase at low bath temperatures, although the system is clearly out of equilibrium. Figure 2(a) indicates that the nonequilibrium phase transition occurs at a temperature which is comparable to the one of the equilibirum model, T eq c = 1.13. However, a closer inspection reveals that the precise value of T c depends on the fixed temperature T 2 in such a way that T c becomes smaller as T 2 is shifted away from the critical value T eq c of the equilibrium model. This conspicuousness is further confirmed by Fig. 2(b), where (for L = 24 and L = 64) m is plotted as function of T for different values of T 2 . Remarkably, the shift of the temperature region (compared to the equilibrium model) where m as function of T increases to large values (indicating the emergence of spin order) is found to be equally large for both system sizes L. This, in turn, signals that the temperature shift of the curves is not a finite size effect (which would vanish for L → ∞), but an actual property of this nonequilibrium vector Potts model. A further interesting observation from Fig. 2(b) is that the nonequilibrium vector Potts model studied here displays an ordered phase, even if one of the heat bath's temperatures is higher than the critical temperature T eq c of the corresponding equilibrium model. However, the critical mean temperature is always below the equilibrium value, as we will discuss below.  Since the spins σ i are coupled to different heat baths, one might expect differences in the phase behavior of the two sublattices L 1 and L 2 . However, as one can see in Fig. 3(b), the Binder cumulant intersects at the same temperature T c in both sublattices. This shows that the transition from the paramagnetic to the ferromagnetic phase occurs collectively in the entire system at the same temperature T c .
For an overview of the critical temperatures in the plane spanned by T 1 and T 2 , we now look at the nonequilibrium phase diagram plotted in Fig. 3(c). The diagonal (black solid) line where T 1 = T 2 corresponds to the equilibrium model. For the nonequilibrium system (T 1 = T 2 ), T c depends on T 1 and T 2 approximately linearly in the vicinity of equilibrium (T 1 ≈ T 2 ) but the dependency becomes strongly nonlinear when T 1 T 2 or T 2 T 1 . This is clearly seen when one compares the actual phase boundary with the dashed curve corresponding to the line along which T = T eq c holds, i.e., T 2 = 2T eq c − T 1 . One can further see that when T 1 = T 2 (i.e., in the equilibrium model), the phase transition occurs at the highest mean temperature T eq c . As soon as there is a temperature difference between the two sublattices, the nonequilibrium phase transitions occur at a lower critical temperature, which deviates from T eq c the more as the difference ∆T = |T 1 − T 2 | between T 1 and T 2 increases. Moreover, there exists a new type of critical temperature T * c = 1.700(2) with the following property: If one sublattice has a temperature higher than T * c , global order is destroyed, irrespective of the temperature of the other sublattice.
Physically, one may understand the phase behavior in the following manner. When heating the system up in the presence of a temperature difference ∆T between the sublattices, disorder is favored already at lower system-averaged temperatures T , showing that a smaller amount of thermal noise destroys the long-range order. Consistent with our physical intuition, a breaking of the translational symmetry (by the temperature difference between the sublattices) reduces the stability of long-range order. Note that this is in sharp contrast to the situation where an homogeneous external magnetic field acts on the system (breaking the up-down symmetry) which increases the stability of long-range ordering.

Critical behavior of specific heat
After the determination of T c and its dependency on the temperatures of the two heat baths, we now turn to the investigation of the thermodynamic properties of our nonequilibrium spin model in the vicinity of the phase transition. First, we calculate the specific heat C v [see Eq. (16)] as function of the mean temperature T for different values of L and T 2 . This quantity is commonly considered in order to characterize second-order phase transitions. As can be seen in Fig. 4(a), C v peaks at a temperature very close to the values of T c that we have previously determined via the Binder cumulant (recall Fig. 3).
For both depicted values of T 2 , the precise location of the peak depends on L in such a way that as the system size is increased, the temperature where the peak is located decreases and approaches T c . We suspect that the peak is exactly at T c in the limit L → ∞, as it is well-known for the equilibrium version of this model. In thermal equilibrium, the specific heat of the model is further known to show universal scaling behavior with respect to the temperature, i.e., with α = 2/3 in the disordered phase [46]. Interestingly, we find that the nonequilibrium model also displays a power-law divergence of C v . Moreover, the critical exponent is the same as in equilibrium, irrespective of the value of the temperature gradient ∆T = |T 2 − T 1 | among the two heat baths (as long as T 1 and T 2 ≤ T * c ). This is exemplarily illustrated for T 2 = 0.3 and T 2 = 0.5 in Fig. 4(b), where C v is plotted for different system sizes (from L = 16 to L = 96) as function of the reduced temperature τ = |1 − T /T c | together with straight (black dashed) lines (with slope −2/3). We checked the scaling behavior for various additional values of T 2 and all of them show a power-law scaling with α = 2/3, demonstrating the robustness of the critical exponent under nonequilibrium conditions. To analyze the critical behavior based on our numerical data in detail, we employ the finite-size scaling technique. To this end, we consider the positions of the peaks of C ν , which give an approximation for the critical temperature as function of the system size L. For the equilibrium 4-state vector Potts model on a square lattice, this quantity scales as ∼ L −ν , with the corresponding critical exponent ν = 2/3 [54]. Also for the nonequilibrium model we obtain ν = 2/3 for all values of T 2 , consistent with the well-known scaling law νd = 2 − α [55] (where d = 2 is the spatial dimension of the lattice).

5.1.2.
Critical behavior of total entropy production Let us now consider the behavior of the total EP which is a direct measure for irreversibility in the sense that it quantifies the distance from equilibrium. To start with, we find that the total EP per spin is always positive, π > 0, whenever T 1 = T 2 . Moreover, π is a convex function of the mean temperature T with minimum at the equilibrium point T = T 1 = T 2 , where π = 0 [consistent with Eq. (7)]. This can be seen in Fig. 5(a), which depicts π vs. T for an exemplary setting of T 2 = 1.5 and L = 32 around the equilibrium mean temperature T = T 1 = T 2 = 1.5.
Depending on whether T 2 is higher or lower than T eq c , the phase transition of the nonequilibrium model lies below or above that minimum (which is always located at T = T 1 = T 2 ). In other words, if T 2 > T c , which is the situation considered in Fig. 5(a), the nonequilibrium phase transition occurs at the left hand side of the minimum, whereas if T 2 < T c , the transition occurs at the right hand side of it. Interestingly, we observe that π as function of T shows a bump around T c = 1.11. In that sense, the function π(T ) itself already signals the occurence of the phase transition.
The dependency of entropy production rate per spin on the temperatures of the two heat baths is plotted in Fig. 5(b), which shows π for different combinations of T 1 and T 2 ranging from 0.1 up to 2.0. As expected, π = 0 whenever T 1 = T 2 [i.e., no temperature gradient ∆T is present and detailed balance is fulfilled, see Eq. (7)]. In contrast to this, there is always a positive rate of entropy production (π > 0) when T 1 = T 2 , consistent with the special case considered in Fig. 5(a). As the gradient ∆T increases, the entropy production rate increases roughly π ∼ ∆T 2 , no matter whether the system is in the PM or the FM phase.
To resolve the EP along individual stochastic trajectories, we plot distributions of the medium entropy production P (φ). Figure 6 displays numerical results for P (φ) at the critical temperature T c [ Fig. 6(a)] and above T c [ Fig. 6(b)], for trajectories of length l = 100 (see Sec. 4 for details). We consider both, the distribution of the entire system as a whole (top panels), and the separate distributions obtained by restricting our observation to one of the two sublattices, L 1 or L 2 , only (middle and bottom panels, respectively). For example, the middle panels show the histograms of all detected values of the medium EP from spins that belong to the sublattice L 1 . Overall, the main characteristics seem to be quite similar for the system at and above the phase transition [compare (a) and (b)]. Let us now take a closer look at the different distributions. Remarkably, in both cases, P (φ) for the whole lattice exhibits a multi-peaked structure [see top panels in Fig. 6]. When inspecting the corresponding sublattice distributions, we notice that the multi-peak structure of the whole system appears to arise as a combination of both sublattices. This is reasonable, as the stochastic trajectories of the whole system expectantly include both, many contributions from the hotter sublattice (which flips more often), and some seldom contributions from the colder sublattice. In fact, the multi-peaked structure looks like a convolution of the distributions from the belonging sublattices. Further, we notice that P (φ) of the individual sublattices have smooth single-peaked shapes. Furthermore, all distributions are discrete, reflecting that the number of possible transitions (and thus, φ values) is finite, because of the discreteness of the underlying spin dynamics. For the colder sublattice, L 2 , φ only takes a particularly small number of values. This is due to the fact that at low bath temperatures, the sublattice only explores a small part of the phase space, and hence, the number of distinct state transitions is small. For the hotter sublattice, L 1 , we notice that the maxima and mean values of P (φ) lie at φ < 0 [in both cases, (a) and (b)]. This alone would violate the second-law of thermodynamics since it implies a negative mean entropy production rate. However, in its usual form, Φ = Π > 0, the second law only applies to the entire system which consists of two sublattices, and the negative mean value simply reflects the heat flows from the hotter to the colder heat bath (overall, the entropy is increased over time).
Next, we study the system size dependency of the total entropy production rate (per spin), π, around the critical point of the phase transition. To this end, we consider a system where T 2 is fixed to a value below T eq c [see the left panel of Fig. 7(a)], and another one where T 2 > T eq c [see the right panel of Fig. 7(a) which is essentially an enlarged version of Fig. 5(a) for different L close to T c ]. Both parts of Fig. 7(a) indicate that π is identical for all system sizes for T values far away from T c , i.e., all lines collapse on a single curves. In contrast, the lines split up around T c , thus, around the phase transition π depends on the system size L. This resembles the behavior of the specific heat [recall Fig. 4]. One further observes the emergence of a shoulder which gets more pronounced while increasing L. It is, however, noteworthy that we do not observe the formation of a saddlepoint or even non-monotonous behavior for all considered system sizes, i.e., until the value L = 96.
In order to study the behavior around the critical point, we inspect the derivative of the entropy production rate dπ/dT for various values of T 2 see Fig. 7. Interestingly, dπ/dT peaks around the temperature of the phase transition. An exception is the case T 2 = T eq c where the total EP naturally vanishes and thus does not peak. Moreover, Figure 7. (a) The entropy production rate per spin, π, as function of the mean temperature T for fixed T 2 and system sizes ranging from L = 16 to L = 96. In the left panel, the temperature of sublattice L 2 is fixed to T 2 = 0.3 which is below the critical temperature T eq c of the equilibrium model, while in the right panel, the temperature of L 2 is T 2 = 1.5, which is above T eq c [see also one observes a dependency of the maximum of dπ/dT on the value of T 2 which (for fixed L) decreases as T 2 approaches T eq c . To analyze the nonequilibrium phase transitions in more detail, we perform a finite-size scaling, similar to our investigation of the specific heat (see Fig. 4). We aim to stress that the application of a finite-size scaling analysis to the EP at a nonequilibrium transition is, to our knowledge, novel. First, we study the scaling behavior of dπ/dT in the disordered phase as function of the reduced temperature τ . Second, we consider the peak height as function of the system size L. As can be seen in Fig. 8(a), dπ/dT shows power-law behavior ∼ τ ζ with an exponent ζ, whose precise value depends on the distance from equilibrium at the phase transition (i.e., on the value of ∆T = |T 2 − T 1 |). Specifically, we detect power-law behavior of dπ/dT for all considered values of T 2 with a decreasing value for ζ as T 2 approaches T eq c , where it nullifies. For T 2 = 0.3 [see the left panel in Fig. 8(a)] the exponent reads ζ = 0.175 (11), while for T 2 = 0.5 [see the right panel in Fig. 8(a)] ζ = 0.145(15) (see the dashed black lines). While the power-law behavior resembles that of the specific heat, there is a marked difference in the sense that the exponent ζ is not constant (such as the exponent α of C v ), but depends on ∆T . In addition, we analyze the scaling behavior of the maximum of dπ/dT as the system size L is increased and show results for T 2 = 0.3 and T 2 = 0.5 in Fig. 8(b). According to the finite-size scaling theory for equilibrium systems [42], all divergent quantities scale as ∼ L a/ν , where a is the critical exponent of the power-law decay of that very quantity. Thus, we test whether the maximum of dπ/dT scales as ∼ L ζ/ν , with ν = 2/3. From our numerical data, we find dπ/dT max ∼ L 0.245 for T 2 = 0.3 and dπ/dT max ∼ L 0.205 for T 2 = 0.5 which is indeed in good agreement with ζ = 0.175(11) (T 2 = 0.3) and ζ = 0.145(15) (T 2 = 0.5) as obtained in Fig. 8(a). The fulfillment of the finite-size scaling relation shows indeed that the derivative of the entropy production rate behaves as a diverging quantity as the critical point of the phase transition is approached. It further demonstrates that the finite-size scaling theory is applicable to the EP rate, despite the dependency of the critical exponent on the temperature gradient between the two sublattices.

BKT-like phase transition in the continuous vector Potts model with q → ∞
Now we turn to the vector Potts model with q → ∞ (also known as the XY model), where the spins can freely rotate in the x − y plane, i.e., all spin orientations σ i ∈ [0, 2π] are allowed. As a consequence of the continuous spin symmetry and the two-dimensional character of the system, there exists no long-range ordered phase at finite temperatures as stated by the Mermin-Wagner theorem [47]. Instead, a quasilong range ordered phase, the BKT phase, occurs at low bath temperatures. While the infinite-order transition between the disordered and the BKT phase is quite well understood in the equilbrium model [56], nonequilibrium BKT phase transitions are in general less understood. In particular, the question of how the EP rate behaves at this transition has, to the best of our knowledge, not been considered in earlier literature. In the previous discussion of the case q = 4, we have seen that the derivative of the total EP shows critical behavior which partially resembles the behavior of the specific heat. Let us now see if this analogy carries over to the BKT transition, which has very different overall characteristics and, in particular, is not accompanied with a divergence of C ν at the critical temperature which is given by T eq c = 0.892880(6) [57] in the equilibrium XY model [56][57][58][59][60][61]. In Fig. 9(a), we show results for π at T 2 = 0.3 and T 2 = 0.5 for system sizes ranging from L = 16 up to L = 64. As indicated there, the EP rate does not split with respect to L in the vicinity of the phase transition. Instead, π is apparently size-independent in the depicted temperature range which includes the BKT transition. In order to visualize the EP rate for different combinations of T 1 and T 2 , we plot π in the T 1 − T 2 plane in Fig. 9(b) together with the derivative of the EP rate with respect to temperature, dπ/dT in Fig. 9(c) for system size L = 32.
Additionally, C v for T 2 = 0.5 and dπ/dT for T 2 = 0.3 and T 2 = 0.5 are plotted in Fig. 10. In contrast to the PM to FM transition of the 4-state vector Potts model, C v in the nonequilibrium XY model does not show any feature like a divergence at criticality. In particular, it only shows a peak around T = 1.1, as does the equilibrium XY model [58], which is above T c . Interestingly, also the derivative of the EP rate with respect to temperature, dπ/dT , does not peak in the vicinity of the critical point. Similar to the specific heat, dπ/dT also shows the peak around T = 1.1 which does not depend on L, i.e., the maximum of dπ/dT does not diverge, but remains constant for all considered system sizes. However, we observe that the maximum of dπ/dT depends on the temperature difference |T 2 − T 1 | between the two sublattices in the vicinity of the peak as confirmed by comparing Fig. 10(b) with Fig. 10(c), where one observes that the maximum value of dπ/dT at T 2 = 0.3 is larger compared to T 2 = 0.5.
Just as for the vector Potts model with q = 4, we investigate the distribution P (φ)  of entropy φ = ∆φ(l) that is produced along stochastic trajectories of length l = 100.
To this end, we plot P (φ) for a system of size L = 64 in the quasi long-range ordered BKT phase at T = 0.5 with T 1 = 0.7 and T 2 = 0.3 (i.e., ∆T = 0.4) in the top panel of Fig. 11(a). The distribution for the whole system seems to be symmetric around the peak position of P (φ) which is located in the positive range, φ > 0 in accordance with the second law of thermodynamics. In contrast, P (φ) for subsystem L 1 peaks in the negative range, and P (φ) for subsystem L 2 peaks at a positive value of φ. This difference in the peak positions just reflects the expected entropy flow from the hot to the cold reservoir. Additionally, one observes different skew directions for P (φ) in the two subsystems. P (φ) for subsystem L 1 is slightly right-skewed, while P (φ) in L 2 is a left-skewed distribution. This effect becomes more pronounced for the system in the PM phase [see Fig. 11(b)] where one clearly observes that P (φ) is skewed in both sublattices. Since the distribution for L 2 is stronger skewed, the distribution for the whole system is also left-skewed.

Conclusions and Outlook
In this paper, we have analyzed the behavior of various critical quantities and that of the total EP rate around the critical point in a nonequilibrium q-state vector Potts model (with q = 4 and q → ∞). The nonequilibrium character results from coupling the spins to two heat baths at different temperatures. Based on this nonequilibrium model, we address several questions: Does the type of phase transition and the critical exponents change by driving the system away from equilibrium? Does the EP exhibit universal behavior around a continuous phase transition? What happens to the EP in the vicinity of an a infinite-order phase transition? First, we have investigated the model with q = 4 in the vicinity of the secondorder phase transition. We found that the critical temperature of the transition decreases as the temperature difference between the two heat baths increases. Moreover, the behavior of the specific heat resembles that of the equilibrium model, i.e., it shows power-law divergence with critical exponents that are independent of the temperature difference. Interestingly, the derivative of the EP rate with respect to temperature behaves, to some extent, similar. It also shows power-law divergence. However, the value of the scaling exponents does depend on the temperature difference and is thus non-universal. Concerning the model with q → ∞, the specific heat as well as the derivative of the EP rate do not show any noticeable behavior around the infinite-order transition from the PM to the quasi long-range ordered BKT phase. Instead, both quantites have a finite peak at a temperature above the critical temperature, i.e., in the PM phase. As the temperature difference between the heat baths increases, the maximum value of the derivative of the EP rate becomes more pronounced. In total, our results provide evidence that the derivative of the EP behaves like a critical quantity, but, as we report here, is non-universal.
Finally, we aim at pointing out perspectives for future work, starting with some questions directly following from the present work. For the sake of generality one should study and compare the behavior of the specific heat with the EP in other dimensions and for different lattice topologies. Further, although the BKT phase transition is not accompanied by a divergence of thermodynamic quantities, in equilibrium it still obeys characteristic scaling dimensions [62]. A more detailed analysis of this transition in the nonequilbrium model, and, specifically, with respect to the derivative of the EP rate, represents an interesting objective of future research.
From a theoretical point of view, it would moreover be worth to think about the connection between EP and specific heat, which seem to behave analogously around criticality, on a fundamental level.
Furthermore, an interesting novel perspective on the nonequilibrium model considered here is the reinterpretation as a model with non-reciprocal coupling between interacing isothermal spins. To be more specific, a vector-Potts model where interacting spins are coupled among each other with two distinct coupling constants (J 1 = J/T 1 and J 2 = J/T 2 ) and uniform temperature follows the exact same equations of motions as our model (with two temperatures and identical coupling constants J). This provides a connection to spin models on directed graphs [63][64][65][66][67][68][69], and to the topic of non-reciprocal interactions, which is currently a focus in nonequilibrium statistical mechanics [70][71][72]. It would be interesting to compare the thermodynamic properties of spin systems subjected to different driving mechanisms, e.g., non-reciprocal couplings, temperature gradients, external fields and colored noise.