Collective power: Minimal model for thermodynamics of nonequilibrium phase transitions

We propose a thermodynamically consistent minimal model to study synchronization which is made of driven and interacting three-state units. This system exhibits at the mean-field level two bifurcations separating three dynamical phases: a single stable fixed point, a stable limit cycle indicative of synchronization, and multiple stable fixed points. These complex emergent dynamical behaviors are understood at the level of the underlying linear Markovian dynamics in terms of metastability, i.e. the appearance of gaps in the upper real part of the spectrum of the Markov generator. Stochastic thermodynamics is used to study the dissipated work across dynamical phases as well as across scales. This dissipated work is found to be reduced by the attractive interactions between the units and to nontrivially depend on the system size. When operating as a work-to-work converter, we find that the maximum power output is achieved far-from-equilibrium in the synchronization regime and that the efficiency at maximum power is surprisingly close to the linear regime prediction. Our work shows the way towards building a thermodynamics of nonequilibrium phase transitions in conjunction to bifurcation theory.


I. INTRODUCTION
While phase transitions are quite well understood at equilibrium, nonequilibrium phase transitions still lack a systematic treatment. They are most commonly described as dynamical phenomenon within the framework of nonlinear dynamics and bifurcation theory [1,2], but their relation to thermodynamics is rarely discussed. This is largely due to the fact that a theory of nonequilibrium thermodynamics was lacking. Stochastic thermodynamics nowadays provides one for systems described by stochastic dynamics [3][4][5][6]. But until recently it has been mostly explored to study noninteracting systems or systems made of few interacting particles. We will use stochastic thermodynamics to explore the physics of nonequilibrium phase transitions in large ensembles of interacting systems.
A motivation to do so which is of great practical importance is to understand how phase transitions, and more generally interactions, affect the performance of large ensembles of nano-machines. Indeed, while these latter have been shown to make very good energy converters, the main drawback remains their low power output. A natural way out is to assemble large numbers of nano-machines, which immediately raises the question of whether certain interactions are favorable to their overall performance. Stochastic thermodynamics provides a powerful framework to do so as it has proved instrumental to analyze the performance of small energy converters operating far-from-equilibrium [3,7,8] (e.g. thermoelectric quantum dots [9,10], photoelectric nanocells [11], molecular motors [12][13][14][15][16]) and their power-efficiency trade-off [17][18][19][20][21]. We emphasize that going beyond linear * Electronic Mail: tim.herpich@uni.lu response is essential here since this is where nonequilibrium phase transitions occur. While some works have been done in this direction, most are restricted to meanfield treatments [22][23][24][25][26][27]. An important aspect of our study will be to analyze in details the emergence of the mean-field description from the underlying stochastic dynamics.
The paradigmatic phase transition which we will consider is synchronization: coupled units with different natural frequencies exhibiting a spontaneous phase-locking to a global frequency [28]. This collective phenomenon was famously described by Huygens [29] who experimentally observed that two pendulum clocks attached to a common support display an "odd kind of sympathy" [29], that is they synchronize in anti-phase. It was later found to be ubiquitous in nature [30]. Synchronization is typically modeled by coupled phase oscillators which exhibit phase-locking when the coupling strength exceeds a critical value [31]. The most commonly used (noisy) Kuramoto model [32][33][34][35] is well understood for an infinite population of oscillator at the mean-field level. Some works also considered few locally coupled oscillators [36][37][38] and even the dissipation resulting from their synchronization [39]. However, little is known about large but finite populations of stochastic oscillators (see e.g. Refs. [40,41]). Progress in this direction was done in Refs. [42][43][44] by introducing a minimal stochastic model made of interacting three-state units and shown to exhibit phase synchronization. It enabled to compare the mean-field dynamics to the Monte-Carlo one. However, since this model is made of three unidirectional stochastic transitions, it does not allow for a consistent thermodynamic description. Furthermore, the extent in which this ingredient is essential for synchronization is not clear. These works also did not provide a detailed understanding of how a linear and irreducible Markov dynamics can give rise to a nonlinear mean-field dynamics with increas-ing system size. This question is particularly intriguing since the Perron-Frobenius theorem ensures that the former dynamics has a unique stationary solution (for finite state spaces) [45] while the latter can exhibit multiple and time-periodic solutions. It is also closely related to the emergence of hydrodynamic modes or metastability [46][47][48][49][50][51][52][53].
In this paper we propose and analyze in great detail a thermodynamically consistent version of the interacting three-state oscillators model. This model can be seen as a toy model for interacting molecular motors [54], enzymes [55,56] or switches [57,58]. At the mean-field level, it displays as a function of the inverse temperature three phases separated by two nonequilibrium phase transitions: a Hopf bifurcation separating a single stable fixed point phase from a stable limit cycle one, and an infinite-period bifurcation separating the limit-cycle phase from a phase with three stable fixed points. At equilibrium only one phase transition survives which separates a phase with a single stable fixed point from one with multiple stable fixed points via a saddle-node bifurcation. A central result is that the spectrum of the Markovian dynamics generator is shown to encode the information about the two bifurcations that are observed in the mean field. The mean-field dynamics is demonstrated to be characterized by the three eigenvalues with dominant real parts (the null one and a complex conjugated pair). At the Hopf bifurcation, a real-part gap between these eigenvalues and the remaining eigenvalues opens up, enabling the emergence of a metastable mean-fieldlike oscillatory dynamics over long times. As the second bifurcation is approached, this difference in real parts further increases while the imaginary parts of the dominant eigenvalues significantly drop causing the oscillations to vanish into three metastable fixed points. The fact that the real part of the most dominant complex conjugated eigenvalue pair converges to zero while the gap with respect to the real parts of all other nonzero eigenvalues becomes larger with increasing system size explains the emergence of the mean-field solutions as the perpetuation of the metastable states. After demonstrating the consistency of stochastic thermodynamics across scales (from the microscopic manybody level to the mean field one), we analyze the dissipated work across the different dynamical regimes. We find that as a function of increasing inverse temperatures the transition towards synchronization is of first order while the outward transition is of second order. A crucial observation is that in the thermodynamic limit, interactions can significantly decrease the dissipated work per oscillator beyond the synchronization threshold and even more so after the second transition towards multistability. Furthermore, interactions in finite assemblies of oscillators enhance this effect in the former case but reduce it in the latter, in particular when the number of oscillators is too low to sustain a long-lasting metastable solution. Finally, we demonstrate that when operating as an energy converter, synchronization significantly enhances the power output per oscillator. Despite operating far-from-equilibrium, the efficiency at maximum power remains quite close to the linear-response prediction of 1/2. Overall, our thermodynamically consistent minimal model for synchronization enables us to reveal with unprecedented detail two complementary facets of a nonequilibrium phase transition: The emergence of different dynamical phases from stochastic dynamics far-from-equilibrium and their thermodynamic characterization using stochastic thermodynamics.
The plan of the paper is as follows. First, in Sec. II, we introduce the description of our model and perform an exact coarse-graining of the dynamics. Next, Sec. III analyzes the different regimes of the mean-field dynamics which motivate the spectral analysis in IV. In Sec. V we compare dynamics between the mean field with finite systems using dynamical Monte Carlo simulations. Furthermore, the thermodynamic laws are formulated in Sec. VI and the work dissipated by noninteracting, small and large interacting networks is compared in VII. Finally, the power-efficiency trade-off in the mean field is investigated in VIII. We conclude with a summary and an outlook to proceeding projects in Sec. IX.

A. Setup
We consider a system consisting of N three-state units with energies i (i = 1, 2, 3). Under the constraint of occupying the same state, units interact globally via an interaction potential. The system is subjected to a nonconservative rotational forcing f and is furthermore in contact with a heat bath at inverse temperature β = (k b T ) −1 , where we set k b ≡ 1 in the following. The schematics of the setup are depicted in Fig. 1.
Owing to the all-to-all interaction, the total interaction energy is obtained by considering a unit in state k and summing up the remaining number of units occupying the same state. It thus holds where u/N is the interaction strength, U 0 = −uN/2 is a constant (∆U 0 = 0) and the notation N k (α) refers to the number of units occupying the single-unit state k in the microstate α. We thus obtain for the change in internal energy The dissipative dynamics of the system, with the above energetics, is described via a Markovian master equation where p α denotes the probability to be in the microstate α. The microscopic transition rates w αα give the probability per unit time for the system to undergo a transition α to α. With only one transition at a time, it follows that the transition rate matrix is irreducible and stochastic, α w αα = 0. This implies the existence of a unique stationary state [45]. We take the microrates to be of Arrhenius form, that is with Γ setting the timescale. The sign function Θ(α, α ) gives preference to transitions down the bias f over counteracting ones. It is defined as Θ(α, α ) = 1 for i (α i − α i ) mod 3 = 1 and Θ(α, α ) = −1 otherwise.
Furthermore, we emphasize that the rates satisfy local detailed balance ensuring that the dynamics is thermodynamically consistent. In the long-time limit t → ∞, the system will tend to its unique steady state, p s α , which is in nonequilibrium due to the presence of the non-conservative driving f .
In absence of driving, microscopic detailed balance holds and along with the local detailed balance condition in Eq. (5) ensures that the equilibrium probability distribution is Gibbsian, i.e. , with the equilibrium free energy Formulating the stochastic process as above gives rise to an exceedingly large state space growing exponentially with the number of units in the network as 3 N − 1. Yet, a closer inspection reveals that a coarse-graining to a mesoscopic space can be done in which the stochastic dynamics can be represented accurately. In fact, the internal energy (and hence also the microscopic transition rates) E(α) ≡ E(N ) does not depend on the topological details encoded in α but only on the mesostate N ≡ (N 1 , N 2 ). The number of microstates α belonging to the same mesostate N is given by if the network is made up of N units. Introducing the marginalized probability P N to observe the mesostate N the ME (3) for the full microstate dynamics can be coarse-grained aṡ with the marginalized mesoscopic transition rates W N N = Ω(N , N )W N N . We note that the coarsegraining preserves the stochastic property and the irreducibility of the transition rate matrix.
The characteristic function χ α ,α emerging in Eq. (11b) is a result of pulling the sums through the microscopic transition rate matrix since the information that w αα = 0, only if transitions between α and α are possible, would be lost. Consequently, the function takes the value 1 if α and α are connected, and is 0 otherwise.
To determine the constrained multiplicity factor Ω(N , N ) ≡ α∈N χ α ,α , we need to address the question of how many microstates α ∈ N are connected with α ∈ N . Two macrostates are connected, if two of three occupation numbers (N 1 , N 2 , N − N 1 − N 2 ) of the macrostates differ by exactly one. In the microstate space, this corresponds to, if, compared entrywise, exactly one number being different in the tuples representing the two microstates. Thus we obtain for the constrained multiplicity factor which indeed does not require any microscopic information. Hence the coarse-graining of the dynamics is exact and leads to a closed ME (11c) represented in terms of mesoscopic states N . This coarse-graining significantly reduces the dimensionality of the state space which grows as [(N + 1)(N + 2)] /2 − 1, thus quadratically as N becomes large. Using Boltzmann's entropy the multiplicity factor of the microstates can be incorporated into the macrorates in a physically meaningful way. The mesoscopic local detailed balance relation thus reads with ∆A(N , N ) = ∆E(N , N )−β −1 ∆S int (N , N ) being the difference in free energy between the macrostates N and N . The mesoscopic sign function Θ(N , N ) is defined analogously to Θ(α, α ). Thus, Θ(N , N ) = 1 for The mesoscopic local detailed balance relation (14) implies that at t → ∞ and for f = 0 the mesoscopic detailed balance holds and the mesoscopic equilibrium probability distribution is again of the Gibbs form with the equilibrium free energy

III. MEAN-FIELD DYNAMICS
In order to further reduce the complexity of the state space of the mesoscopic ME (11c) we first operate in the mean-field (MF) limit where N → ∞. In this limit, the total change in internal energy due to a transition in Eq. (2b) simplifies and the corresponding scaled current density J(n i , n j ) ≡ lim N →∞ W N N /N becomes where n i = N i /N denotes the occupation density of the single-unit state i and Θ(i, j) = 1 for (i − j) mod 3 = 1, while Θ(i, j) = −1 otherwise. The evolution equation for the mean occupation density reads In the MF approximation we replace any n-point correlation function with a product of n averages thus yieldinġ which represents a closed nonlinear equation. The validity of this approximation can be proved [45] in the macroscopic limit N → ∞. Hence the MF system can be described by a single three-state unit, where the (average) occupation density of the single-unit states is assigned to the three states of the MF unit. We therefore identify the MF occupation density, n i , as the probability for any unit to occupy the single-unit state i = 1, 2, 3. Its dynamics is ruled by the nonlinear MF equatioṅ with the MF transition rates obeying local detailed balance Unit conservation erases one degree of freedom such that there are only two independent variables n 1 and n 2 . We proceed by choosing a flat energy landscape, i.e. by setting i = const ∀i. This allows us to immediately read off the symmetric point n * i = 1/3 as an analytic solution to the nonlinear MF Eq. (21). Linearizing the Eq. around this fixed point (FP),ṅ i = j ∂kij ∂nj nj =n * j n j , we find for the eigenvalues of the linear stability matrix For attractive interactions (u < 0) between the units the real part of λ ± changes its sign at β c1 = −3/u. This crossover suggests that the stable symmetric FP destabilizes and degenerates into a limit cycle (LC) corresponding to a Hopf bifurcation indicative of synchronization.
In appendix A, the LC is characterized in the vicinity of the Hopf bifurcation which is shown to occur supercritical for attractive interactions. Moreover, a closer inspection of the MF rates in Eq. (23) reveals that for any f and β repulsive interactions, u > 0, always lead to the stable symmetric FP. Fig. 2 depicts the MF phase space for different β and f in units of u. The symmetric fixed point is only stable for β < β c1 . We find in agreement with Eq. (24) that for finite f there is a phase characterized by stable LCs if β ≥ β c1 . For any value of f , there is an additional phase with three non-symmetric FPs for β ≥ β c2 (f ). We set u = −1 in the following and briefly address a subtlety of the MF system. In Fig. 2 the analytic solution to Eq. (21), n i = 1/3, is chosen as initial condition. In fact, at temperatures close to the first critical temperature β c1 the long-time solution is initial-condition dependent: For 0 < f < f c ≈ 0.21 there is a finite set of initial conditions different from the symmetric FP that will not lead to a LC but to a non-symmetric stable fixed point. If f ≥ f c , the dynamics will always exhibit a LC regardless of the chosen initial condition.
Before studying the different nonequilibrium phases of this model, we discuss it for f = 0, i.e. at equilibrium. Figure 3a) shows, starting from the initial condition n(0) = (1/3, 1/3) , the long-time solution n eq 1 (t) for different values of β.
At the critical temperature β c1 the system exhibits three non-symmetric stable FPs that emerge via a saddlenode bifurcation. Our thermodynamic framework allows us also to work within the nomenclature of statistical mechanics. Interestingly, the saddle-node bifurcation corresponds to a first-order equilibrium phase transition since the derivative of the MF free energy with respect to β at the critical point β c1 is divergent. Starting from the symmetric FP, these attractive FPs are observed to move towards the corners of the triangle in the n eq 1 −n eq 2 plane. This is physically plausible since at low temperatures the system tends to occupy its lowest energy state where all units are occupying the same state. The dependence of the multiple equilibrium states on the initial condition in the low-temperature phase is investigated in Fig. 3b).
In the lower triangle, n eq 1 is plotted as a function of all physical initial conditions (n 1 (0),n 2 (0)). As a complement, the other MF probability n eq 2 is shown in the upper triangle, where the axis labels are omitted for better readability. Each triangle exhibits two phases which are separated by a contour line. Combining these two panels [59], we find that for every physical initial condition the system will eventually arrive at one of the three non-symmetric FPs. These differ only by permutations of their components (n eq 1 , n eq 2 , 1 − n eq 1 − n eq 2 ), where two of them are identical according to the two phases in each of the panels in b). Figure 4a) depicts the MF probability n 1 (t) as a function of β at long times for f = 1.0 at which the range of β values for which LCs can be observed is close to its maximum, according to Fig. 2. In agreement with Eq. (24), the oscillations emerge at β c1 for any finite f . The oscillations exhibit an increasing frequency with β up to a point where they slow down. At the second critical point, β c2 (f = 1.0) ≈ 6.11, the oscillation period diverges corresponding to an infinite-period bifurcation [60]. The initial-condition dependence of the stationary states in the non-symmetric asynchronous phase (NA) is depicted in d), with β = 7.0. Again, depending on the chosen initial condition, the system will eventually arrive in one of the three non-symmetric FPs, which are again related to each other by permutations of their components. Here, in contrast to the equilibrium case, all components are different. This reflects the presence of the force distorting the symmetry of the states. The distortion occurs since it is more likely to jump from the largely populate state into the lower occupied state following the bias rather than the opposite way. This asymmetry naturally increases as the system is driven further out-of-equilibrium. This constitutes the first important result. We have developed a minimal model which, according to Eqs. (24) and (A16), exhibits synchronization and is thermodynamically consistent due to Eqs. (5), (14) and (23). We also note that synchronization only occurs in a finite range of temperatures: Fig. 3 illustrates that at low temperatures the equilibrated system is energydriven and tends to its energetic ground state, while for very high temperatures the system is entropy-driven and takes a uniform stationary probability distribution. By extrapolation from equilibrium to the non-equilibrium scenario where the synchronization phase emerges, we realize that Fig. 4 invites for an analogous physical interpretation of the low-and high-temperature limit in the non-equilibrium case. Moreover, the limit β → 0 represents equilibrium since forward and backward transition for each pair of states becomes equally probable for any f and thus detailed balance holds. We remark furthermore that the term "minimal" refers to the dimensionality of the MF dynamics given by Eq. (21), which is a natural requirement to observe synchronization since a single-variable nonlinear differential equation can naturally not have complex eigenvalues.

IV. SPECTRAL ANALYSIS: METASTABILITY
A crucial aspect of our model is that it allows us to study its (thermo-)dynamic features for large but finite system sizes and in particular to monitor the convergence of the stochastic dynamics to the MF dynamics. In order to proceed, we begin by stating the formal solution to the mesoscopic ME (11c) that reads where P (0) is the initial probability distribution, λ i are the eigenvalues and Φ L i , Φ R i are the left-and right eigenvectors of the non-symmetric real transition rate matrix W constituting an orthonormal dual basis The index i * characterizes, if existent, the modes with eigenvalues and eigenvectors being the complex-conjugated to those labeled with i. The Perron-Frobenius theorem (PFT) [45] stipulates that for this irreducible, autonomous and stochastic matrix there is a non-degenerate eigenvalue, the Perron-Frobenius eigenvalue (PFE), λ 0 = 0, which is strictly greater than the real part of any other eigenvalue, Re(λ i ) < λ 0 ∀i = 0. Note that the labeling of the eigenvalues is given by the order of their real parts, 0 > Re(λ 1 ) > . . . . Consequently, Eq. (25b) has a unique, infinite-time solution, P s = c 0 Φ R 0 , characterized by the PFE and the associated right eigenvector Φ R 0 . Hence the stationary state of the mesoscopic system P s cannot exhibit stable oscillations (S phase) or multistability (NA phase). On the other hand, one expects that the transition from the mesoscopic system to the MF is smooth as the system size N grows. This apparent paradox is caused by the non-commutation of the infinitetime limit t → ∞ and the mean-field limit N → ∞, i.e.
The right-hand side corresponds to the symmetric stationary state of the SA phase for all temperatures, while the left-hand side is temperature dependent: For β c1 ≤ β < β c2 the system is in a time-periodic state (S phase) and for β ≥ β c2 the dynamics will go to one of the nonsymmetric steady states (NA phase) depending on the chosen initial condition. At β < β c1 the left-hand side also corresponds to the symmetric stationary state, hence the two limits commute only at sufficiently high temperatures.
To resolve this apparent contradiction we look for clues in the spectrum of the Markov generator in the mesoscopic ME (11c) and establish a link between finite-size systems and MF via the notion of metastability. Even though the stationary state is inevitably reached in the infinite-time limit, there could be long-living metastable states that display the phenomenology of the MF. The time-scales to characterize such a state are encoded in the spectrum as follows where τ r is the relaxation time to reach the unique steady state, i.e. it specifies the time-scale at which all finitetime modes have been removed from the dynamics. τ m is the metastable time at which all modes have decayed except for those forming the metastable state, that is the one associated with the eigenvalue λ 1 and the stationary one characterized by the PFE λ 0 . Here, we assume that only a pair of modes associated with a complex-conjugated non-null eigenvalue is contributing to metastability, while there could be an arbitrary number of modes forming the metastable state. This assumption will be numerically verified in the following. Physically, τ l corresponds to the lifetime of that metastable state. To reconcile the stochastic dynamics with its asymptotic solution in the macroscopic limit, the MF dynamics, τ l is required to become increasingly larger with the system size N , while τ m remains finite since the different MF phases emerge at finite time. Using Eqs. (27a)-(27c), these prerequisites translate into conditions on the real parts of the dominant eigenvalues of the Markov generator: The real-part gap between the two first non-null eigenvalues, Re(λ 1 ) − Re(λ 2 ), has to increase by Re(λ 1 ) converging to zero (corresponding to a diverging relaxation time τ r ), while Re(λ 2 ) has to approach a finite value (assuring the emergence of the metastable phenomena at finite times). Moreover to mimic MF dynamics the metastable state has to be oscillatory (Im(λ 1 ) = 0) in the S phase and quasistationary (Im(λ 1 ) = 0) in the NA phase.
Before addressing the question of how the stochastic dynamics converges to the MF, we study the real parts a) and the imaginary parts b) of the two dominant nonzero eigenvalues of the spectrum in all three different phases (2 ≤ β ≤ 8) for a system size of N = 300 in Fig. 5. We remark that for all β, these two eigenvalues in fact occur as complex-conjugated pairs and only those with a positive imaginary part are depicted in panel b).
Furthermore, to stress that the different phases of the finite-size system for β > β c1 are only present for finite times, we rename them differently than in the MF: asynchronous phase (A), synchronous metastable phase (SM) and asynchronous metastable phase (AM).  Figure 5. The real part a) and the imaginary part b), as well as the ratio, of the two most dominant eigenvalues with distinct real part, λ1,λ2, with positive imaginary part are depicted as a function of β. In addition, the LC frequency ω lc that is numerically extracted from the asymptotic (t → ∞) MF dynamics is compared to the imaginary part of the most dominant eigenvalue. All eigenvalues correspond to a generator W for a system of size N = 300. Panel c) shows the lifetime of the metastable state τ l as function of β for different system sizes. The labels of the different phases, that is the asynchronous phase (A), the synchronous metastable phase (SM) and the asynchronous metastable phase (AM) are in correspondence with the labels of the different phases in the MF limit introduced in the preceding Sec. III.
As can be seen in panel a), the real parts of the two eigenvalues both approach zero up to β ≈ 4 followed by a monotonic decrease of Re(λ 2 ) while Re(λ 1 ) changes only slightly and for β > β c2 rapidly goes to zero. According to Eq. (27c), this observation along with the fact that Re(λ 1 )/ Re(λ 2 ) drops at both critical points (dashed lines) suggests that the lifetime τ l of the metastable state is increasing from the SM to the AM regime. The small values of | Re(λ 1 )| in the SM and AM phase and the sharp changes in the ratio of the real parts at both critical points provide a first hint that the metastable state is constituted by only the stationary mode and those associated with the first complex-conjugated non-null eigenvalue.
This claim is further strengthened by studying the corresponding imaginary parts of these eigenvalues as shown in Fig. 5b). We find an excellent agreement in the SM phase between the LC frequency ω lc in the MF that is numerically extracted from the dynamics and Im(λ 1 ). The LC frequency ω lc only coincides with the imaginary part of the Jacobian from the linear stability analysis in Eq. (24) at the bifurcation point β c1 , where the linearization of the nonlinear ME (21) is exact. Moreover, the ratio between the imaginary parts of λ 1 and λ 2 remains nearly constant at 0.5 within the A and SM phase implying that the frequency of oscillation of the mode corresponding to λ 2 is half as that of λ 1 . In the AM phase Im(λ 1 ) quickly goes to zero consistent with our MF observations that show no oscillations.
Consistent with the discussion of the real parts, Fig.  5c) illustrates that the lifetime of the metastable state is nearly zero in the A phase and starts to increase significantly at the first critical point up to a local maximum in the SM phase. The lifetime τ l is monotonically decreasing for larger β before it sharply rises in the AM phase. All clues thus indicate that in the two phases where the MF exhibits non-unique solutions at infinite times, the finite system displays metastability. As expected, for all temperatures in the metastable phases the lifetime is monotonically increasing with N .
Next, to shed some light on the convergence of the finite-system dynamics to the MF dynamics, we investigate the changes in the spectrum as we approach the MF limit. To this end, we look at the first few dominant non-zero eigenvalues as a function of the system size N at β = 4 representing the SM phase. We observe in Fig. 6a) that the real parts of these eigenvalues are approaching the PFE. Though the inset reveals an increasing time-scale separation between the mode associated with λ 1 and the faster decaying modes for larger systems. The monotonically increasing behavior of τ l and τ r with N implies an increasing lifetime of the metastable state, while this time window is shifted to increasingly larger times, hence the finite-system dynamics are converging to the MF. To be fully consistent with the MF, the metastable state must be appearing in the dynamics at a finite time. Taking into account all the aforementioned hints (encoded in Fig. 5 and to be made in the following) that indeed only the modes associated with λ 1,1 * are contributing to the metastability and therefore in correspondence with the MF solution, it is reasonable to expect that Re(λ 2 ) converges to a finite value for larger N . Unfortunately, extracting the dominant eigenvalues of the generator for even larger N is not feasible.
As another striking evidence for the hypothesis that the metastable state comprises only the stationary and the first non-null complex-conjugated mode, the imaginary part of the dominant eigenvalue λ 1 smoothly converges to the LC frequency ω lc in the MF while the imaginary parts of other modes display a distinct separation as seen in Fig. 6b). This is confirmed in Figs. 6c) -f) depicting the mean occupation densities, n(t) = N N /N P N (t), using the full spectral decomposition of the Markov generator in Eq. (25b) and the truncated one for N = 10 2 , 10 3 at β = 4.
To understand the metastability in the AM phase, Fig.  7 depicts the real and imaginary parts of the eigenvalues associated with the most dominant modes in panels a) and b), respectively, as a function of N for β = 7. In contrast to Fig. 6a), here, Re(λ 2 ) clearly converges to a finite value with Re(λ 1 ) quickly going to zero already for small N . This is confirmed by the inset showing that τ l and τ r take very large values already for smaller systems implying that the metastability in the AM phase is much stronger than in the SM phase. As expected, in compliance with the nonoscillatory MF solution, the small magnitudes of the imaginary part vanish rapidly with growing system size as displayed in panel 7b). Figs. 7c) -d) reaffirm that the metastable state in the AM phase is reached at short time-scales and is quasistationary. Moreover, we note the large time-scales (cf. the scale of the axis of the insets) over which the metastable state can be observed in the dynamics in compliance with the observations made in panel 7a). Thus, we conclude from the observations made in this section, that for sufficiently large systems in the SM and AM phase at times τ m t τ r , the relaxation dynamics is determined by the metastable state associated with λ 1,1 * and the PFE. This time span is increasing with N [cf. Figs 6a) and 7a)] such that the metastable states can be observed over increasingly larger times. Owing to the PFT, any finite system will eventually leave these metastable states at times t τ r and relax into the unique stationary state at infinite time. To sum up, we obtain the important result that the different phases and bifurcations of the MF dynamics are encoded in the spectrum of the Markov generator.

V. SIMULATIONS
Solving the ME (11c) for systems on the order of N ∼ 10 3 via full diagonalization of the propagator is computationally not feasible [61]. Hence for extremely large systems we resort to a stochastic simulation algorithm for computing the time evolution of the (Markov) jump processes. This dynamic Monte Carlo method, often referred to as Gillespie algorithm [62,63], generates trajectories of a stochastic process that are exact solutions to the stochastic process. By generating sufficiently many trajectories one can infer the statistics of the observables of the stochastic process, in particular the average values generically denoted by · . Figure 8 depicts the n 2 − n 1 plots generated with the Gillespie algorithm sampling over 10 6 trajectories for selected values of β and for different system sizes N = 10 2 , 10 4 . Except for β = 6.1 shown in e), the larger system, N = 10 4 , agrees well with the MF limit at the displayed times. The smaller system, N = 10 2 , significantly deviates in both the SM phase (β = 4, 5, 6.1) and AM phase (β = 7). In the A phase (β = 2, 3), there are no visible differences between the different finite system sizes and the MF limit solution, as all are relaxing into the unique symmetric fixed point [red closed circle in panel a)]. Of particular interest is the dynamics for β = 7. While the smaller system directly goes to the stationary state, the larger system quickly approaches and wiggles around the FP of the MF limit. This can be seen from the inset that displays a magnification around one of the MF FP [orange closed circle in f)]. Depending on the initial condition the metastable state will approach one of the three MF FPs.
This shows that the stochastic dynamics of sufficiently large systems indeed reproduces the MF dynamics at long times and thus confirms all predictions made above based on the spectral analysis. As an exception, we observe in e) that close to the infinite-period bifurcation, β ≈ β c2 , the large system does not exhibit the characteristics of the solution in the MF limit. However, an even larger system, N = 10 6 , shows signatures of the LC albeit still deviating. These deviations are due to the strong fluctuations in the vicinity of the phase transition calling for larger N such that the finite system can accurately represent the deterministic limit. We remark that this feature is also manifested in the increasing deviations between the LC frequency, ω lc , and the imaginary part of the crucial eigenvalue, λ 1 , as the second critical point, β c2 ≈ 6.1068 is approached [cf. Fig. 5c)].
However, there is a set of initial conditions for which the stochastic dynamics will not go to one of these metastable states. This set of initial conditions is readily constructed via all possible linear combinations of right eigenvectors of the mesoscopic generator from Eq. (11c), P (0) = i =1 a i Φ R i , excluding the mode associated with the crucial eigenvalue pair λ 1,1 * . It follows from the orthonormal dual-basis property of the eigensystem that the weights c 1,1 * = 0 in Eq. (25b). Hence the metasta- bility would be removed from the dynamics.
This prompts the question whether the metastability is a generic (up to a negligibly small set of special initial conditions) property of the stochastic process or just an artifact of choosing suitable initial conditions. This question is addressed in Fig. 9, where the initial conditions are sampled and the joint probability distribution P {n 1 (t = 20), n 2 (t = 20)} for different system sizes N = 10 2 , 10 4 and β = 4, 7 is shown in a density plot. In panel 9a) the distribution exhibits its maxima indicated by the red spots close to the corners of the LC in the MF limit. Overall, the distribution clearly exhibits signatures of the LC but the probability mass is still dispersed around the LC contour. Moreover, over the entire state space there are regions with finite probability. If the system size is notably increased to N = 10 4 , as depicted in Fig. 9b), the probability mass is sharply concentrated on the LC contour.
Turning to panels c) and d) corresponding to the AM regime with β = 7, we observe that the joint probability distribution for the smaller system already reproduces to a good approximation the three non-symmetric FPs in the MF limit. The distribution for the larger system further concentrates the probability mass on the three FPs as can be seen by comparing the insets on the left and on the right magnifying the vicinity of the FPs. The convergence of the probability distribution at smaller N to the MF limit for larger β is consistent with the observations already made in the spectral analysis in Fig.  5. We thus confirm, once again, that the metastability and therefore the convergence to the MF limit increases with N and β. Next, and more importantly, the emergence of the metastable state(s) is, up to a negligible set of special initial conditions, indeed a generic property of the stochastic process. It is insightful to monitor the time evolution of P {n 1 (t), n 2 (t)} starting from a uniform grid at t = 0 up to a time as the distribution becomes stationary or time-periodic. To this end, the supplementaries [64] include movies displaying the dynamics of the distributions shown in Fig. 9.
We have so far established a connection between linear stochastic dynamics and the deterministic nonlinear MF dynamics via the study of the spectrum of the Markov generator. Indeed, the different dynamical phases and bifurcations in the MF are encoded in the spectrum and appear as metastable states for long times in the stochastic dynamics. These predictions are confirmed by our simulations. We now proceed by analyzing the bifurcations as nonequilibrium phase transitions in the thermodynamic observables. In doing so, we link deterministic bifurcation theory to stochastic thermodynamics.

VI. THERMODYNAMIC LAWS
We first introduce the basic thermodynamic state functions in this model: the microscopic internal energy and the system entropy For our setup with an autonomous driving, f , these functions can only change due to the time-dependence of the probability distribution. The rate of change of internal energy naturally defines the microscopic first law of thermodynamics with the heat and work rate given by where the sign function Θ(α, α ) is defined below Eq. (4). The microscopic local detailed balance relation (5) can be expressed in terms of the heat exchange with the bath along the forward transition The system entropy change d t s = ṡ e + σ (34) can be decomposed into the entropy flow from the bath to the system and the non-negative entropy production (EP) rate Equation (36) is the second law of thermodynamics and the inequality follows straightforwardly from ln x ≤ x−1.
The marginalization of the microscopic probability p α performed in Sec. II, yet being exact on the level of the dynamics, does not a priori guarantee that the thermodynamic observables defined above are invariant under this coarse-graining [65]. Defining E N to be the internal energy of the system in the macrostate N , and applying the coarse-graining from Eq. (10) on the expression for the internal energy in Eq. (29a), we obtain Thus the coarse-graining admits a representation in the mesospace while it keeps the internal energy invariant. The heat and work fluxes can also be exactly coarsegrained as Consequently, the first law of thermodynamics has a closed mesoscopic representation which is identical to the one from Eq. (30). We note that after the coarse-graining the heat increment is no longer directly given by the local detailed balance relation like in the microspace, cf. Eq. (33), but also contains the internal entropy from Eq. (13) [65]. We define the system entropy in the mesospace consisting of the non-equilibrium entropy defined by Eq. (29b) and the internal entropy accounting for the multiplicity of distinct microscopic configurations for a given macrostate. Analogously to Eq. (34), we decompose the time-derivative of the entropy into the entropy flow and the EP rate The definitions in Eqs. (40), (42) are in general not coinciding with those made at the microscopic level, i.e. S = s , Σ = σ . The nonlinearity of the system entropy and the EP [Eqs. (29b), (36)] in the microstate probability p α is incompatible with the coarse-graining. Instead, an application of Eq. (36) gives rise to additional entropic contributions which are dependent on microscopic information, hence the coarse-grained equation can not be closed [65]. For the special case of a stationary probability distribution, P s N , one can show (cf. appendix B) via the spanning tree formula [66] that the microstates belonging to the same macrostate are equally probable, p α = P N /Ω(N ). In the stationary limit, the entropies in mesoscopic representation are therefore identical to those in microscopic representation, i.e. S s = s s , Σ s = σ s . For this particular case, the second law boils down to the steady entropy flow Ṡ s e being equal to the magnitude of the steady EP rate Σ s . Using the non-positivity of the average stationary heat, β Q s ≤ 0, we easily verify that Σ s ≥ 0.
We now turn to the MF case and consistently define the first law in this limit with the heat and work fluẋ where i, j = 1, 2, 3 specifies the state of the single MF unit. In analogy to Eq. (29b), we write the system entropy in the MF limit as which we split into the MF entropy floẇ and the non-negative MF EP ratė As the MF represents the asymptotic limit of the mesospace, it holds that all the mesoscopic averages of the intensive observables X /N that are consistent with the coarse-graining in Eq. (10) converge to the corresponding observables X in the MF limit, lim N →∞ Ẋ N =Ẋ , with X = E, Q, W, S e , N i . Consequently, for the MF definitions in Eqs. (46) and (48) to represent the physical entropies, we have to restrict to the stationary case, n s , which yields for the second law in the MF limiṫ The non-negativity of the MF EP follows from the nonpositivity of the MF heat in this model. We have thus developed three different levels (microspace, mesospace and MF) to consistently characterize the energetics of our model. For the first law, the lower levels of description are equivalent, while for the second law they only coincide in the stationary limit. The same applies asymptotically in the macroscopic limit to the thermodynamic observables defined at the MF level.

VII. DISSIPATED WORK
With the thermodynamic framework developed in the preceding section at hand, we can now proceed by addressing one of the crucial research questions of this work, that is the thermodynamics of non-equilibrium phase transitions. We are naturally interested in the (metastable) synchronization regime bounded by the two phase transitions. Since the nonstationary EP represented in the microscopace is not identical to the one in the mesospace [Eqs. (36) and (42)], we characterize the nonequilibrium phase transitions via the dissipated work given by Eqs. (38b) and (45b). At metastable or infinite time, the work is observed to be always dissipative on average, that is the system takes up the energy from the nonconservative force, W > 0, and dissipates it into the bath in the form of heat, Q < 0, for all temperatures and system sizes. Figure 10a) depicts the difference between the stationary work current of a single unit, W 1 = 2 Γf sinh(f β/2), and the asymptotic work current per unit in a network of size N , W N ≡ W N t as a function of β for different N . The derivation of the single-unit stationary work current, W 1 , is deferred to appendix C and given by Eq. (C6). The asymptotic work current, W N , is numerically determined by solving Eqs. (11c) and (21) for a finite and a MF system, respectively. As seen in Fig. 10a), the large (N = 10 4 ) system agrees excellently with the MF limit for all temperatures, while the smaller systems, albeit showing a qualitatively similar behavior, unlike the dynamics, deviate significantly.
Since the single unit work current is governed by a smooth and convex function, we observe that the dissipated MF work exhibits striking changes at the critical points β c1,2 . The vicinities of these critical points are magnified in the two insets. The phase transitions in the dissipated MF work at β c1 and β c2 exhibit a kink and a saddle, respectively, and are therefore reminiscent of a first-and second-order equilibrium phase transition. Remarkably, owing to the metastability in the stochastic dynamics, sufficiently large systems also exhibit finitetime signatures of these nonequilibrium phase transitions at the bifurcation points which blur out with decreasing system size.
In the high-temperature limit, β → 0, the difference ∆W 1N ≡ W 1 − W N between the dissipated work of a single unit and an interacting system per unit is always zero since the interaction energy gets canceled (N i = N j in Eq. (2b) and u/N → 0 as N → ∞.). While for the MF this holds true in the entire SA phase, for finite systems the range of β values in the A phase for which the interaction energy is negligible decreases with N . We find that interactions reduce the costs to maintain the system in its nonequilibrium state, ∆W 1N > 0. This work dissipation gap, ∆W 1N , is a monotonically increasing function of β and becomes infinitely large in the lowtemperature limit, since ∆W 1N /W 1 → 1 as β → ∞. This asymptotic limit can be seen as follows. In Sec. III we observed that in the low-temperature limit, one can make use of the equilibrium picture where the system tends to occupy its energy ground states. In this limit, we have for the dissipated work of a finite network per unit Hence we have shown that at low and intermediate temperatures an interacting network of any size is energetically favorable with respect to a noninteracting one. Interestingly, in the the two phases of higher temperature, the operational costs per unit can be further de- creased by employing smaller networks. As one approaches the second critical point the different curves intersect and in the NA/AM phase the operation of larger networks gives rise to less work dissipation per unit. This is also illustrated in Fig. 10b) that depicts the difference in the dissipated work between a system of size N = 10 4 exhibiting metastability and a smaller sys-tem which does not display metastable states. In agreement with panel a), the smaller system requires less input per unit to be maintained in the two higher temperature phases, since the difference ∆W N,10 4 ≡ W N − W 10 4 < 0 while the opposite holds true in the AM phase, where ∆W N,10 4 > 0.
Again, we observe at the critical points significant changes in ∆W N,10 4 : At the first critical point ∆W N,10 4 takes a local minimum and at the second critical point it changes sharply around an inflection point. It is plausible that these changes are more pronounced for decreasing N as the reference system (N = 10 4 ) exhibits metastability, such that for increasing differences in the network size compared the distance to metastable behavior implying signatures of phase transitions in the dissipated work becomes larger.
For the same reasons as stated in the context of plot 10a), ∆W N,10 4 goes to zero in the high-temperature regime, while in the low-temperature limit one obtains if N < 10 4 . This limit is illustrated by the purple closed circle in the plot. For the larger system the work difference is decreasing in the range of available data. Generating data for larger β to monitor the convergence to the low-temperature limit is not possible since the simulation becomes numerically unstable owing to the large values the exponentials take in the transition rates.
To illustrate the data underlying the plots in Fig. 10b), we show in panel c) to e) the time-scaled work asymptotics per unit for different system sizes as chosen for the blue curve in panel b) as well as the MF limit for selected values of β = 2, 4, 7. We note the excellent agreement between the MF limit and the large system in compliance with the observations made in panel 10a). On the other hand, the small system clearly deviates from the large systems in all three different regimes, even though we observed that in the SA/A phase the dynamics of large and small systems can hardly be distinguished. Due to the approximate time-periodicity in the S/SM phase, the dissipated work is also oscillating.
Finally, Fig. 11 depicts the difference between the stationary single-unit and the asymptotic MF unit work current, ∆W 1∞ , as a function of β for different f . Again, ∆W 1∞ = 0 in the A phase since the single and the MF unit are indistinguishable in the high-temperature regime as shown above in the context of Fig. 10a). For β ≥ β c1 the second critical point is gradually shifting to smaller β [cf. Fig. 2] while the difference ∆W 1∞ is monotonically increasing with decreasing f . Therefore, if compared to the MF, the additional costs to maintain the nonequilibrium stationary state of the noninteracting system at a given temperature are the smaller the further it is driven out-of-equilibrium. This implies in particular that the dissipation of the synchronized system at fixed temperature is approaching the one of the non-synchronized system as they are further driven out-of-equilibrium. To summarize, we have obtained two major results in this section. First, though the nonequilibrium phase transitions are naturally only present in the MF-limit where the nonlinear dynamics exhibits the supercritical Hopf and the infinite-period bifurcation, we find that the metastability observed in the finite-system dynamics translates into signatures of genuine nonequilibrium phase transition. This consistently connects linear stochastic dynamics, nonlinear deterministic dynamics, and thermodynamics and furthermore demonstrates that thermodynamics of nonequilibrium phase transitions and bifurcation theory are closely related. Secondly, any finite and attractive interaction in a network reduces the dissipated work per unit. Interestingly, if operating in the synchronous phase, it is even more economic to employ interacting but smaller networks. What is still open to investigate is how the nonequilibrium phase transitions affect the power-efficiency trade-off, if the system operates as an energy-converting machine.

VIII. EFFICIENCY AT MAXIMUM POWER
In order to construct such an energy converter with our system both a positive force f 1 > 0 and a negative force f 2 < 0 are applied on the same unit. Examples for this type of work-to-work conversion are could be double quantum dot channel capacitively coupled to a quantum point contact [10] or the biological motors kinesin and myosin. In the latter case, the motor is driven forward with f 1 by extracting energy via ATP hydrolysis while the load carried by the motor is modeled as f 2 [26,67]. In general, these two forces obey two different distributions accounting for the crucial fluctuations these motors exhibit. Since the following discussion is restricted to the MF limit, we consider the homogenous case where the same positive and negative force are applied on all units.
We thus decompose the net force f = f 1 + f 2 into the driving force f 1 > 0 and the load force f 2 < 0. Their respective steady-state work contributions are denoted by W s 1 and W s 2 . Substituting Eq. (23) into Eq. (49), yields the following decomposition of the stationary EP in the MF limit where S W s k i = β W s k , k = 1, 2. Based on Eq. (53), we use as an unambigious definition for the efficiency of this work-to-work conversion (cf. Refs. [68,69]) At equilibrium (f = 0), the reversible limit, η c = 1 is attained while out of equilibrium (f = 0) the efficiency is bounded, 0 < η < 1. Of particular interest is the efficiency at maximum power (EMP) [70], which results from the optimization of the stationary output power P ≡ ∂W s 2 /∂t with respect to the output force The maximization parameter f * is determined by the condition ∂P/∂f 2 = 0, while fixing f 1 = 1 and thus varying the total dissipation. In the SA phase, β < β c , the stationary power putput coincides with the average work current of a single unit given by Eq. (C6). For the other two phases (S and NA), we have to resort to simulations to obtain the power output. Moreover, owing to the time-periodic state in the S phase, the power is periodically changing in time. Hence we consider the time-average of the power over one LC period. Figure 12a) shows the numerically determined output power P as a function of β and f 2 in a density plot.
The white dashed lines indicate the critical points β c1,2 as a function of the output force. Thus the area enclosed by those lines corresponds to the S phase. Remarkably, we find that the maximum output power is generated in this phase. In particular, the global maximum of the output power indicated by the purple closed circle lies inside the S phase. At large β that represents the NA phase, the generated power rapidly drops. In panel b) the output power maximized with respect to the output force for different values of the inverse temperature is depicted. The numerical data from panel a) is overlaid with the (semi-)analytic results in the SA phase (green solid line) and the low-temperature limit and shows an excellent agreement. These limiting cases can be obtained as follows. In the SA phase, the condition for maximization of the power results in a transcendental equation that must be treated numerically. In the low-temperature limit, the extremum Figure 12. Depiction of the output power a) as a function of the output force f2 and the inverse temperature β. The white dashed lines correspond to the numerically determined critical points as a function of the output force. Hence the enclosed area defines the synchronization phase S. The global maximum of the output power is indicated by the purple closed circle. In panel b) the maximum output power P * is optimized with respect to f2 and in panel c) the associated EMP η * (f * 2 ) is displayed. In panel c) the dashed lines specify the critical points and the synchronization phase S. The efficiency at the global maximum power is indicated by the purple arrow. The (semi-)analytic solution for β < βc [green lines] is overlaid with the numerical data in the lower panels.
can not be satisfied for any f 2 compatible with the constraint β = ∞. The efficiencies associated with the processes corresponding to the data points in panel 12b) are depicted in panel c). Again, the semianalytic solution for the temperatures corresponding to the SA phase (green solid line) is compared with the numerical results and shows an excellent agreement at these temperatures. As β approaches zero, the EMP takes the universal linearresponse value for tightly-coupled (only one net-current) systems, η * = 0.5 η c [71,72]. This can be seen by expanding the expression for the stationary work current in the SA phase given by Eq. (C6) up to first order in β which yields the linear-response relation J s ≈ L f with the Onsager coefficient L = Γ β. Therefore, small products βf correspond to linear response in our model and lead to EMP values very close to 1/2. With increasing β, the system starts to respond nonlinearly and the efficiency decreases monotonically and nonlinearly.
It is worth emphasizing that the efficiency for the global maximum power output achieved in the far-fromequilibrium S phase and indicated by the purple closed circle is still close to the universal linear-response EMP value. This finding points out the importance of nonequilibrium phase transitions for the performance of an assembly of nano-machines and suggests synchronization as an operating mode faciliating very efficient energyconversion processes with appreciable power output.

IX. CONCLUSION AND PERSPECTIVES
We introduced and studied a thermodynamically consistent minimal model of N driven and globally interacting three-state units obeying linear Markovian dynamics.
The mean-field dynamics (which is exact when N → ∞) exhibits two nonequilibrium phase transitions as a function of the inverse temperature, a Hopf and an infinite-period bifurcation. These separate three distinct phases consisting respectively of a stable fixed point where all units states are equiprobable, a limit cycle corresponding to synchronization of the units, a coexistence of three stable fixed points where the units states have unequal probabilities.
We demonstrated that these transitions are encoded in the spectrum of the generator of the linear Markovian dynamics. The two dominant complex-conjugated eigenvalues, beside the null one, describe the mean-field dynamics over metastable times (i.e. times located between the inverse of the real parts of the next dominant eigenvalues and the inverse of their own real part) which increase with N . All predictions based on the spectral analysis were confirmed employing dynamic Monte Carlo simulations.
After having established a nonequilibrium thermodynamics description of our model at different scales, we characterized the nonequilibrium phase transitions using the work dissipated by the external force driving the units. The mean-field dissipated work which reproduces very well the large N results undergoes a first order phase transition followed by a second order one as a function of the inverse temperature. When comparing a single unit to a unit in an interacting network, the average dissipated work for both units is equal in the first phase, while for the interacting unit it remarkably drops in the synchronization phase and drops even further in the third phase. Interestingly, in the presence of interactions and when N is too low to produce a meaningful metastable mean-field dynamics, the average dissipated work in the second (resp. third) phase is lower (resp. higher) than for in the mean field (N → ∞).
Finally, when operating our system in the mean-field limit as a work to work converter, we found that the synchronization phase leads to a significant boost in the power output. The efficiency at maximum power of this far-from-equilibrium machine is surprisingly close to the universal linear-regime prediction.
The model we used is minimal in that it contains the minimal ingredients to be thermodynamically consistent and at the same time give rise to a limit cycle. As most minimal stochastic models, it may find various applications (e.g. interacting molecular motors or coupled quantum dots). The methods we used are generic in that they can be used on other models.
A natural extension of this work would consist in analyzing thermodynamic fluctuations in particular close to phase transitions based on generating function techniques and large deviation theory. Another one would be to explore the effects of local interactions and of the network topology on the dissipated work. While the qualitative behavior of synchronization is likely to survive [42,43], new interesting spatiotemporal regimes may emerge [73].
At the fundamental level, our work shows an instance where bifurcation theory can be augmented with a thermodynamic interpretation to move towards a theory of nonequilibrium phase transitions. In such a theory, bifurcations would arise from the nonlinearities of the mean-field dynamics which emerges from an underlying stochastic thermodynamics of interacting systems in the macroscopic limit.
From a more utilitarian perspective, our work suggest interesting avenue towards engineering interactions between assemblies of small machines to efficiently generate power, in particular in far-from-equilibrium regimes where nonequilibrium phase transitions may arise.
where z is a complex variable, z * its complex-conjugate, ∆β = β − β c1 gives the distance of the inverse temperature to the critical inverse temperature of the Hopf bifurcation and g = O z 2 is a smooth function of (z, z * , ∆β).
Appendix B: Equal-probability of stationary microstates belonging to a macrostate A special case for which also the EP and system entropy can be exactly represented by macrostate ensemble quantities is the nonequilibrium steady state reached at large times. The probabilities associated with states in the stationary regime can be calculated via the spanning tree formula. We denote the graph representing the network by G. A spanning tree, T (G) of a graph is defined as a covering subgraph of G, i.e. all of its edges are also edges of G and it contains all vertices (microstates) of G. It is furthermore required that T (G) is connected and contains no circuits. We introduce the notation A(T (µ) α (G)) referring to the µth spanning tree rooted in α, that is a tree whose branches are pointing towards the vertex α. The spanning tree formula states [66] p s α = µ A(T As was already discussed above, the transition rates do not depend on the microstates belonging to the same pair of macrostate. Moreover, the connectivity of the network is also not a function of the microstate, since, due to the all-to-all interaction, the number of edges of any vertex in the microspace network is always 2N , such that the number of spanning trees rooted in α is constant for all α inside the same macrostate. Thus, at steady state, all microstates constituting the same macrostate where Ω(N ) is the number of microstates forming the macrostate N given by a trinomial coefficient of the occupation numbers N i determined in Eq. (9).

Appendix C: Stationary solution for single unit
We consider a single unit with states i = 1, 2, 3 whose evolution is governed by the ME where P is the (macro-)probability to find the unit in the single state i with the transition rates with the sign function Θ(i, j) as defined in Eq. (23) ensuring the validity of local detailed balance. The steadystate work current reads Using the spanning tree formula from Eq. (B1), one obtains for the stationary probabilities P s 1 = a 1 a 1 + a 2 + a 3 , P s 2 = a 2 a 1 + a 2 + a 3 , For a flat energy landscape, i =const, we indeed find that the symmetric stationary solution P i = 1/3 is independent of β and f like in the MF limit. Next, the stationary work current is given by