# Collective power:
Minimal model for thermodynamics of nonequilibrium phase transitions

Tim Herpich  Electronic Mail<EMAIL_ADDRESS>Juzar Thingna  Massimiliano
Esposito  Complex Systems and Statistical Mechanics, Physics and Materials
Science Research Unit, University of Luxembourg, L-1511 Luxembourg, Luxembourg

(August 27, 2024)

###### Abstract

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 Strogatz (2001); Guckenheimer and
Holmes (1983), 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 Seifert (2012); den Broeck and Esposito (2015); Zhang
_et al._ (2012); Rao and Esposito (2018). 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 Seifert
(2012); Pekola (2015); Ciliberto (2017) (e.g. thermoelectric quantum dots
Esposito  _et al._ (2009a); Cuetara and Esposito (2015), photoelectric
nanocells Rutten  _et al._ (2009), molecular motors Ge  _et al._ (2012);
Gaspard and Gerritsma (2007); Lau  _et al._ (2007); Seifert (2011); Altaner
_et al._ (2015)) and their power-efficiency trade-off Shiraishi  _et al._
(2016); Proesmans  _et al._ (2016); Pietzonka and Seifert (2018); Lee and Park
(2017); Polettini and Esposito (2017). We emphasize that going beyond linear
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 mean-field treatments Vroylandt  _et al._ (2017); Shiraishi
(2017); Golubeva and Imparato (2012, 2012, 2013); Imparato (2015); Sasa
(2015). 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 Pikovsky  _et al._ (2003). This collective
phenomenon was famously described by Huygens Huygens (1895) who experimentally
observed that two pendulum clocks attached to a common support display an “odd
kind of sympathy” Huygens (1895), that is they synchronize in anti-phase. It
was later found to be ubiquitous in nature Strogatz (2004). Synchronization is
typically modeled by coupled phase oscillators which exhibit phase-locking
when the coupling strength exceeds a critical value Winfree (1967). The most
commonly used (noisy) Kuramoto model Kuramoto (2003); Strogatz (2000); Acebrón
_et al._ (2005); Gupta  _et al._ (2014) is well understood for an infinite
population of oscillator at the mean-field level. Some works also considered
few locally coupled oscillators Sakaguchi  _et al._ (1987); Daido (1988);
Walgraef  _et al._ (1983) and even the dissipation resulting from their
synchronization Izumida  _et al._ (2016). However, little is known about large
but finite populations of stochastic oscillators (see e.g. Refs. Risler  _et
al._ (2004, 2005)). Progress in this direction was done in Refs. Wood  _et
al._ (2006a, b, 2007) 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 increasing 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) Van Kampen (2007) while the
latter can exhibit multiple and time-periodic solutions. It is also closely
related to the emergence of hydrodynamic modes or metastability Gaspard (1996,
1998); Gaveau and Schulman (1998); Biroli and Kurchan (2001); Esposito and
Gaspard (2005a, b); Macieszczak  _et al._ (2016); Rose  _et al._ (2016).

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 Campàs  _et
al._ (2006), enzymes Hill (1977); Hill and Chen (1981) or switches Raymo  _et
al._ (2003); Simão  _et al._ (2011). 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-field-like 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.

##  II Model

###  II.1 Setup

We consider a system consisting of $N$ three-state units with energies
$\epsilon_{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 non-conservative rotational forcing $f$ and is furthermore in contact
with a heat bath at inverse temperature $\beta=(k_{b}T)^{-1}$, where we set
$k_{b}\equiv 1$ in the following. The schematics of the setup are depicted in
Fig. 1.

Figure 1: Small network of identical and all-to-all interacting three-state
units in contact with a heat bath $\beta$ and in the presence of a
nonconservative force $f$.

We denote a microstate by the multiindex
$\alpha$=$(\alpha_{1},\ldots,\alpha_{i},\ldots,\alpha_{N})$ with
$\alpha_{i}\\!=\\!1,2,3$. As an example, ordering the units from left to right
and from top to bottom, the microstate displayed in Fig. 1 reads
$\alpha$=$(2,1,3,2,3,1,2,2,3)$. Consider a transition from microstate
$\alpha^{\prime}$ to $\alpha$ amounting to a transition between the single
unit energy states $\epsilon_{j}$ to $\epsilon_{i}$. For such a transition the
occupation numbers change according to $N_{i}\\!\to\\!N_{i}+1$ and
$N_{j}\\!\to\\!N_{j}-1$. To determine the change in internal energy $\Delta
E(\alpha,\alpha^{\prime})\\!=\\!\epsilon(\alpha)\\!-\epsilon(\alpha^{\prime})\\!+U(\alpha)\\!-U(\alpha^{\prime})$,
the total interaction energy $U(\alpha)$ of the network is required. 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

$\displaystyle
U(\alpha)=\frac{u}{N}\sum\limits_{k=1}^{3}\;\sum\limits_{l=1}^{N_{k}\\!(\alpha)\\!-\\!1}\\!l=\frac{u}{2N}\sum\limits_{k=1}^{3}N_{k}^{2}(\alpha)+U_{0},$
(1)

where $u/N$ is the interaction strength, $U_{0}\\!=\\!-uN/2$ is a constant
$(\Delta U_{0}\\!=\\!0)$ and the notation $N_{k}(\alpha)$ refers to the number
of units occupying the single-unit state $k$ in the microstate $\alpha$. We
thus obtain for the change in internal energy

$\displaystyle\Delta E(\alpha,\alpha^{\prime})\\!$
$\displaystyle=\\!\epsilon(\alpha)\\!-\\!\epsilon(\alpha^{\prime})\\!+\\!\frac{u}{2N}\\!\sum\limits_{k=1}^{3}\\!\left[N_{k}^{2}(\alpha)\\!-\\!N_{k}^{2}(\alpha^{\prime})\right]$
(2a) $\displaystyle=\epsilon_{i}-\epsilon_{j}+\frac{u}{N}(N_{i}-N_{j}+1).$
(2b)

###  II.2 Master-Equation

The dissipative dynamics of the system, with the above energetics, is
described via a Markovian master equation (ME)

$\displaystyle\dot{p}_{\alpha}=\sum\limits_{\alpha^{\prime}}w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}\,,$
(3)

where $p_{\alpha}$ denotes the probability to be in the microstate $\alpha$.
The microscopic transition rates $w_{\alpha\alpha^{\prime}}$ give the
probability per unit time for the system to undergo a transition
$\alpha^{\prime}$ to $\alpha$. With only one transition at a time, it follows
that the transition rate matrix is irreducible and stochastic,
$\sum_{\alpha}w_{\alpha\alpha^{\prime}}=0$. This implies the existence of a
unique stationary state Van Kampen (2007). We take the microrates to be of
Arrhenius form, that is

$\displaystyle w_{\alpha\alpha^{\prime}}$
$\displaystyle=\Gamma\;\mathrm{e}^{-\frac{\beta}{2}\left(\Delta
E(\alpha,\alpha^{\prime})-\Theta(\alpha,\alpha^{\prime})f\right)}\,,$ (4)

with $\Gamma$ setting the timescale. The sign function
$\Theta(\alpha,\alpha^{\prime})$ gives preference to transitions down the bias
$f$ over counteracting ones. It is defined as
$\Theta(\alpha,\alpha^{\prime})\\!=\\!1$ for
$\sum_{i}(\alpha_{i}\\!-\\!\alpha^{\prime}_{i})\text{ mod }3=\\!1$ and
$\Theta(\alpha,\alpha^{\prime})\\!=\\!-1$ otherwise. Furthermore, we emphasize
that the rates satisfy local detailed balance

$\displaystyle\ln\frac{w_{\alpha\alpha^{\prime}}}{w_{\alpha^{\prime}\alpha}}=-\beta\left(\Delta
E(\alpha,\alpha^{\prime})-\Theta(\alpha,\alpha^{\prime})f\right)\,,$ (5)

ensuring that the dynamics is thermodynamically consistent. In the long-time
limit $t\to\infty$, the system will tend to its unique steady state,
$p^{s}_{\alpha}$, which is in nonequilibrium due to the presence of the non-
conservative driving $f$.

In absence of driving, microscopic detailed balance

$\displaystyle
w_{\alpha\alpha^{\prime}}\,p^{eq}_{\alpha^{\prime}}=w_{\alpha^{\prime}\alpha}\,p^{eq}_{\alpha}\,,$
(6)

holds and along with the local detailed balance condition in Eq. (5) ensures
that the equilibrium probability distribution is Gibbsian, _i.e._ ,

$\displaystyle p^{eq}_{\alpha}$
$\displaystyle=\mathrm{e}^{-\beta\left(E_{\alpha}-A^{eq}\right)}\,,$ (7) with
the equilibrium free energy $\displaystyle A^{eq}$
$\displaystyle=-\frac{1}{\beta}\,\ln\,\sum\limits_{\alpha}\mathrm{e}^{-\beta
E_{\alpha}}\,.$ (8)

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(\alpha)\\!\equiv\\!E(\bm{N})$ does not depend
on the topological details encoded in $\alpha$ but only on the mesostate
$\bm{N}\\!\equiv\\!(N_{1},N_{2})$. The number of microstates $\alpha$
belonging to the same mesostate $\bm{N}$ is given by

$\displaystyle\Omega(\bm{N})\equiv\sum\limits_{\alpha\in\bm{N}}1=\binom{N}{N_{1}}\binom{N-N_{1}}{N_{2}}=\frac{N!}{\prod\limits_{i}N_{i}!}\,,$
(9)

if the network is made up of $N$ units. Introducing the marginalized
probability $P_{\bm{N}}$ to observe the mesostate $\bm{N}$

$\displaystyle P_{\bm{N}}$
$\displaystyle\equiv\sum\limits_{\alpha\in\bm{N}}\,p_{\alpha}\,,$ (10)

the ME (3) for the full microstate dynamics can be coarse-grained as

$\displaystyle\dot{P}_{\bm{N}}$
$\displaystyle=\sum\limits_{\alpha\in\bm{N}}\sum\limits_{\bm{N}^{\prime}}\sum\limits_{\alpha^{\prime}\in\bm{N}^{\prime}}w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}$
(11a)
$\displaystyle=\sum\limits_{\bm{N}^{\prime}}\tilde{W}_{\bm{N}\bm{N}^{\prime}}\sum\limits_{\alpha\in\bm{N}}\sum\limits_{\alpha^{\prime}\in\bm{N}^{\prime}}\;p_{\alpha^{\prime}}\,\chi_{\alpha^{\prime},\alpha}$
(11b)
$\displaystyle=\sum\limits_{\bm{N}^{\prime}}W_{\bm{N}\bm{N}^{\prime}}\;P_{\bm{N}^{\prime}}\,,$
(11c)

with the marginalized mesoscopic transition rates
$W_{\bm{N}\bm{N}^{\prime}}=\Omega(\bm{N},\bm{N}^{\prime})\tilde{W}_{\bm{N}\bm{N}^{\prime}}$.
We note that the coarse-graining preserves the stochastic property and the
irreducibility of the transition rate matrix.

The characteristic function $\chi_{\alpha^{\prime},\alpha}$ emerging in Eq.
(11b) is a result of pulling the sums through the microscopic transition rate
matrix since the information that $w_{\alpha\alpha^{\prime}}\neq 0$, only if
transitions between $\alpha^{\prime}$ and $\alpha$ are possible, would be
lost. Consequently, the function takes the value 1 if $\alpha^{\prime}$ and
$\alpha$ are connected, and is 0 otherwise.

To determine the constrained multiplicity factor
$\Omega(\bm{N},\bm{N}^{\prime})\equiv\sum_{\alpha\in\bm{N}}\,\chi_{\alpha^{\prime},\alpha}$,
we need to address the question of how many microstates
$\alpha\\!\in\\!\bm{N}$ are connected with
$\alpha^{\prime}\\!\in\\!\bm{N}^{\prime}$. 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

$\displaystyle\Omega(\bm{N},\bm{N}^{\prime})$
$\displaystyle=N^{\prime}_{1}\,\delta_{N^{\prime}_{1},N_{1}+1}+N^{\prime}_{2}\,\delta_{N^{\prime}_{2},N_{2}+1}+$
$\displaystyle+\\!(N\\!\\!-\\!\\!N^{\prime}_{1}\\!\\!-\\!\\!N^{\prime}_{2})\,\delta_{N\\!-(N^{\prime}_{1}+N^{\prime}_{2}),N\\!-(N_{1}\\!+\\!N_{2})+1}\,,$
(12)

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 $\bm{N}$. This coarse-graining significantly
reduces the dimensionality of the state space which grows as
$\left[(N+1)(N+2)\right]/2-1$, thus quadratically as $N$ becomes large. Using
Boltzmann’s entropy

$\displaystyle S^{int}(\bm{N})=\ln\Omega(\bm{N})\,,$ (13)

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

$\displaystyle\ln\frac{W_{\bm{N}\bm{N}^{\prime}}}{W_{\bm{N}^{\prime}\bm{N}}}=-\beta\left[\Delta
A(\bm{N},\bm{N}^{\prime})-\Theta(\bm{N},\bm{N}^{\prime})f\right],$ (14)

with $\Delta A(\bm{N},\bm{N}^{\prime})\\!=\\!\Delta
E(\bm{N},\bm{N}^{\prime})\\!-\\!\beta^{-1}\Delta
S^{int}(\bm{N},\bm{N}^{\prime})$ being the difference in free energy between
the macrostates $\bm{N}^{\prime}$ and $\bm{N}$. The mesoscopic sign function
$\Theta(\bm{N},\bm{N}^{\prime})$ is defined analogously to
$\Theta(\alpha,\alpha^{\prime})$. Thus, $\Theta(\bm{N},\bm{N}^{\prime})=1$ for
$\sum_{i}(\bm{N}_{i}\\!-\bm{N}_{i}^{\prime})\text{ mod }3=\\!1$ and
$\Theta(\bm{N},\bm{N}^{\prime})\\!=\\!-1$ otherwise.

The mesoscopic local detailed balance relation (14) implies that at
$t\to\infty$ and for $f=0$ the mesoscopic detailed balance

$\displaystyle
W_{\bm{N}\bm{N}^{\prime}}\,P^{eq}_{\bm{N}^{\prime}}=W_{\bm{N}^{\prime}\bm{N}}\,P^{eq}_{\bm{N}}\,,$
(15)

holds and the mesoscopic equilibrium probability distribution

$\displaystyle
P^{eq}_{\bm{N}}=\mathrm{e}^{-\beta\left(A_{\bm{N}}-A^{eq}\right)},$ (16)

is again of the Gibbs form with the equilibrium free energy

$\displaystyle
A^{eq}=-\frac{1}{\beta}\ln\sum\limits_{\bm{N}}\mathrm{e}^{-\beta
A_{\bm{N}}}\,.$ (17)

##  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\\!\to\\!\infty$. 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})\\!\equiv\\!\lim_{N\to\infty}W_{\bm{N}\bm{N}^{\prime}}/N$
becomes

$\displaystyle
J(n_{i},n_{j})=\Gamma\,\mathrm{e}^{-\frac{\beta}{2}\\!\left(\epsilon_{i}-\epsilon_{j}+u(n_{i}-n_{j})-\Theta(i,j)f\right)}\,n_{j}\,,$
(18)

where $n_{i}\\!=\\!N_{i}/N$ denotes the occupation density of the single-unit
state $i$ and $\Theta(i,j)=1$ for $(i-j)\text{ mod }3=\\!1$, while
$\Theta(i,j)\\!=\\!-1$ otherwise. The evolution equation for the mean
occupation density reads

$\displaystyle\langle\dot{n}_{i}\rangle$ $\displaystyle=\sum\limits_{j\neq
i}\langle J(n_{i},n_{j})\rangle-\langle J(n_{j},n_{i})\rangle.$ (19)

In the MF approximation we replace any $n$-point correlation function with a
product of $n$ averages thus yielding

$\displaystyle\dot{\overline{n}}_{i}\equiv\langle\dot{n}_{i}\rangle$
$\displaystyle=\sum\limits_{j\neq i}J(\langle n_{i}\rangle,\langle
n_{j}\rangle)-J(\langle n_{j}\rangle,\langle n_{i}\rangle),$ (20)

which represents a closed nonlinear equation. The validity of this
approximation can be proved Van Kampen (2007) in the macroscopic limit
$N\\!\leavevmode\nobreak\ \to\leavevmode\nobreak\ \\!\infty$. 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,
$\overline{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 equation

$\displaystyle\dot{\overline{n}}_{i}\\!$
$\displaystyle=\sum\limits_{j}k_{ij}\,\overline{n}_{j}\,,$ (21)

with the MF transition rates

$\displaystyle
k_{ij}=\Gamma\,\mathrm{e}^{-\frac{\beta}{2}\\!\left(\epsilon_{i}-\epsilon_{j}+u(\overline{n}_{i}-\overline{n}_{j})-\Theta(i,j)f\right)}\,,$
(22)

obeying local detailed balance

$\displaystyle\ln\frac{k_{ij}}{k_{ij}}$
$\displaystyle=-\beta\left(\epsilon_{i}-\epsilon_{j}+u(\overline{n}_{i}-\overline{n}_{j})-\Theta(i,j)\,f\right).$
(23)

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
$\epsilon_{i}\\!=\\!\operatorname{const}\,\forall i$. This allows us to
immediately read off the symmetric point $\overline{n}_{i}^{*}=1/3$ as an
analytic solution to the nonlinear MF Eq. (21). Linearizing the Eq. around
this fixed point (FP),
$\dot{\overline{n}}_{i}=\sum\limits_{j}\left.\tfrac{\partial
k_{ij}}{\partial\overline{n}_{j}}\right|_{\overline{n}_{j}=\overline{n}_{j}^{*}}\overline{n}_{j}$,
we find for the eigenvalues of the linear stability matrix

$\displaystyle\lambda_{\pm}\\!=\\!-\Gamma(\beta
u\\!+\\!3)\cosh\left(\frac{\beta
f}{2}\right)\\!\pm\\!\mathrm{i}\sqrt{3}\,\Gamma\sinh\left(\frac{\beta
f}{2}\right).$ (24)

For attractive interactions ($u<0$) between the units the real part of
$\lambda_{\pm}$ changes its sign at $\beta_{c_{1}}=-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
$\beta$ repulsive interactions, $u>0$, always lead to the stable symmetric FP.

Fig. 2 depicts the MF phase space for different $\beta$ and $f$ in units of
$u$. The symmetric fixed point is only stable for $\beta<\beta_{c_{1}}$. We
find in agreement with Eq. (24) that for finite $f$ there is a phase
characterized by stable LCs if $\beta\geq\beta_{c_{1}}$. For any value of $f$,
there is an additional phase with three non-symmetric FPs for
$\beta\geq\beta_{c_{2}}(f)$.

Figure 2: Phase space in the MF varying the parameters $\beta$ and $f$ in
units of $u$. The black lines correspond to the set of critical points
$\beta_{c_{1}}$ and $\beta_{c_{2}}(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), $\overline{n}_{i}=1/3$,
is chosen as initial condition. In fact, at temperatures close to the first
critical temperature $\beta_{c_{1}}$ the long-time solution is initial-
condition dependent: For $0<f<f_{c}\approx 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\geq 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 $\bm{\overline{n}}(0)=(1/3,1/3)^{\top}$, the long-time
solution $\overline{n}^{eq}_{1}(t)$ for different values of $\beta$.

Figure 3: Density plot of the equilibrium occupation density
$\overline{n}^{eq}_{1}$ for different $\beta$ and times $t$ for an initial
condition $\bm{\overline{n}}(0)$ equal to the symmetric FP in a), and as a
function of all physical initial conditions $\overline{n}_{1}(0)$ and
$\overline{n}_{2}(0)$ at time $t=10^{3}$ and for $\beta=4.0$ in b) (lower left
triangle). For completeness, the upper right triangle in panel b) shows the
other component $\overline{n}^{eq}_{2}$. The times are chosen to be
sufficiently long such that the system has relaxed to equilibrium.

At the critical temperature $\beta_{c_{1}}$ the system exhibits three non-
symmetric stable FPs that emerge via a saddle-node 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 $\beta$ at the critical point $\beta_{c_{1}}$ is
divergent. Starting from the symmetric FP, these attractive FPs are observed
to move towards the corners of the triangle in the
$\overline{n}^{eq}_{1}\\!-\\!\overline{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, $\overline{n}^{eq}_{1}$ is plotted as a function of all
physical initial conditions ($\overline{n}_{1}(0)$,$\overline{n}_{2}(0)$). As
a complement, the other MF probability $\overline{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 111Note that a state does not correspond to a folding of the
two triangles but a rotation of one of the two planes about $180^{\circ}$ and
subsequent overlapping of the two layers., 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
$(\overline{n}^{eq}_{1},\overline{n}^{eq}_{2},1-\overline{n}^{eq}_{1}-\overline{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 $\overline{n}_{1}(t)$ as a function of
$\beta$ at long times for $f=1.0$ at which the range of $\beta$ 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 $\beta_{c_{1}}$ for any
finite $f$. The oscillations exhibit an increasing frequency with $\beta$ up
to a point where they slow down. At the second critical point,
$\beta_{c_{2}}(f=1.0)\approx 6.11$, the oscillation period diverges
corresponding to an infinite-period bifurcation Keener (1981). The initial-
condition dependence of the stationary states in the non-symmetric
asynchronous phase (NA) is depicted in d), with $\beta=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.

Figure 4:  Illustration of the occupation probability $\overline{n}_{1}$ as a
function of $\beta$ and $t$ for $f\\!=\\!1.0$ in a), while b) shows the
occupation densities $\overline{n}_{1}$ in the lower left triangle and
$\overline{n}_{2}$ in the upper right triangle as a function of all initial
conditions at $\overline{n}_{1}(0)$ and $\overline{n}_{2}(0)$ at time
$t\\!=\\!10^{3}$ and for $\beta\\!=\\!7.0$. The initial condition underlying
the plots in panels a) is $(\overline{n}_{1}(0)=1,\overline{n}_{2}(0)=0)$. The
times are chosen such that the system has reached either the unique FP in the
symmetric asynchronous phase (SA), the stable LC in the synchronous phase (S),
or one of the three non-symmetric FPs in the non-symmetric asynchronous phase
(NA). The triple points defined by the intersecting contour lines in b)
correspond to the symmetric unstable FP present in the NA phase.

This constitutes the first important result. We have developed a minimal model
which, according to Eqs. (24) and (73), 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 energy-driven
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 $\beta\to
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

$\displaystyle\bm{P}(t)$ $\displaystyle=\mathrm{e}^{\bm{W}t}\,\bm{P}(0)$ (25a)
$\displaystyle=\\!\sum\limits_{i,i^{*}}\\!\mathrm{e}^{\lambda_{i}t}\\!\underbrace{\Phi^{L}_{i}\bm{P}(0)}_{\equiv
c_{i}}\Phi^{R}_{i}+c_{i^{*}}\mathrm{e}^{\lambda_{i^{*}}t}\Phi^{R}_{i^{*}},$
(25b)

where $\bm{P}(0)$ is the initial probability distribution, $\lambda_{i}$ are
the eigenvalues and $\Phi^{L}_{i}$, $\Phi^{R}_{i}$ are the left- and right
eigenvectors of the non-symmetric real transition rate matrix $\bm{W}$
constituting an orthonormal dual basis $\Phi^{L}_{i}\Phi^{R}_{j}=\delta_{ij}$.
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) Van Kampen (2007) stipulates that for this
irreducible, autonomous and stochastic matrix there is a non-degenerate
eigenvalue, the Perron-Frobenius eigenvalue (PFE), $\lambda_{0}=0$, which is
strictly greater than the real part of any other eigenvalue,
$\operatorname{Re}(\lambda_{i})<\lambda_{0}\;\forall i\neq 0$. Note that the
labeling of the eigenvalues is given by the order of their real parts,
$0>\operatorname{Re}(\lambda_{1})>\ldots\;$. Consequently, Eq. (25b) has a
unique, infinite-time solution, $\bm{P}^{s}=c_{0}\Phi^{R}_{0}$, characterized
by the PFE and the associated right eigenvector $\Phi_{0}^{R}$.

Hence the stationary state of the mesoscopic system $\bm{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 infinite-time limit $t\to\infty$ and the mean-field
limit $N\to\infty$, _i.e._

$\displaystyle\lim_{t\to\infty}\lim_{N\to\infty}\;\bm{P}(t)\neq\lim_{N\to\infty}\;\underbrace{\lim_{t\to\infty}\;\bm{P}(t)}_{\bm{P}^{s}},\quad\text{if
}\;\beta\geq\beta_{c_{1}}.$ (26)

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 $\beta_{c_{1}}\leq\beta<\beta_{c_{2}}$ the system is in a time-periodic
state (S phase) and for $\beta\geq\beta_{c_{2}}$ the dynamics will go to one
of the non-symmetric steady states (NA phase) depending on the chosen initial
condition. At $\beta<\beta_{c_{1}}$ 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

$\displaystyle\tau_{r}$
$\displaystyle\sim-\frac{1}{\operatorname{Re}(\lambda_{1})}$ (27a)
$\displaystyle\tau_{m}$
$\displaystyle\sim-\frac{1}{\operatorname{Re}(\lambda_{2})}$ (27b)
$\displaystyle\tau_{l}$
$\displaystyle\equiv\tau_{r}-\tau_{m}\sim\frac{1}{|\operatorname{Re}(\lambda_{1})|}\left(1-\frac{\operatorname{Re}(\lambda_{1})}{\operatorname{Re}(\lambda_{2})}\right),$
(27c)

where $\tau_{r}$ is the relaxation time to reach the unique steady state,
_i.e._ it specifies the time-scale at which all finite-time modes have been
removed from the dynamics. $\tau_{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 $\lambda_{1}$ and the stationary one
characterized by the PFE $\lambda_{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, $\tau_{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, $\tau_{l}$ is required to become
increasingly larger with the system size $N$, while $\tau_{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,
$\operatorname{Re}(\lambda_{1})-\operatorname{Re}(\lambda_{2})$, has to
increase by $\operatorname{Re}(\lambda_{1})$ converging to zero (corresponding
to a diverging relaxation time $\tau_{r}$), while
$\operatorname{Re}(\lambda_{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
($\operatorname{Im}(\lambda_{1})\neq 0$) in the S phase and quasistationary
($\operatorname{Im}(\lambda_{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
non-zero eigenvalues of the spectrum in all three different phases
($2\leq\beta\leq 8$) for a system size of $N=300$ in Fig. 5. We remark that
for all $\beta$, 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
$\beta>\beta_{c_{1}}$ 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,
$\lambda_{1}$,$\lambda_{2}$, with positive imaginary part are depicted as a
function of $\beta$. In addition, the LC frequency $\omega_{lc}$ that is
numerically extracted from the asymptotic ($t\to\infty$) MF dynamics is
compared to the imaginary part of the most dominant eigenvalue. All
eigenvalues correspond to a generator $\bm{W}$ for a system of size
$N=\leavevmode\nobreak\ 300$. Panel c) shows the lifetime of the metastable
state $\tau_{l}$ as function of $\beta$ 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 $\beta\approx 4$ followed by a monotonic decrease of
$\operatorname{Re}(\lambda_{2})$ while $\operatorname{Re}(\lambda_{1})$
changes only slightly and for $\beta>\beta_{c_{2}}$ rapidly goes to zero.
According to Eq. (27c), this observation along with the fact that
$\operatorname{Re}(\lambda_{1})/\operatorname{Re}(\lambda_{2})$ drops at both
critical points (dashed lines) suggests that the lifetime $\tau_{l}$ of the
metastable state is increasing from the SM to the AM regime. The small values
of $|\operatorname{Re}(\lambda_{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 $\omega_{lc}$ in the MF
that is numerically extracted from the dynamics and
$\operatorname{Im}(\lambda_{1})$. The LC frequency $\omega_{lc}$ only
coincides with the imaginary part of the Jacobian from the linear stability
analysis in Eq. (24) at the bifurcation point $\beta_{c_{1}}$, where the
linearization of the nonlinear ME (21) is exact. Moreover, the ratio between
the imaginary parts of $\lambda_{1}$ and $\lambda_{2}$ remains nearly constant
at $0.5$ within the A and SM phase implying that the frequency of oscillation
of the mode corresponding to $\lambda_{2}$ is half as that of $\lambda_{1}$.
In the AM phase $\operatorname{Im}(\lambda_{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 $\tau_{l}$ is monotonically decreasing for larger
$\beta$ 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 $\beta=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 $\lambda_{1}$ and the faster
decaying modes for larger systems. The monotonically increasing behavior of
$\tau_{l}$ and $\tau_{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 $\lambda_{1,1^{*}}$ are contributing to the metastability and therefore
in correspondence with the MF solution, it is reasonable to expect that
$\operatorname{Re}(\lambda_{2})$ converges to a finite value for larger $N$.
Unfortunately, extracting the dominant eigenvalues of the generator for even
larger $N$ is not feasible.

Figure 6: Real a) and corresponding imaginary b) parts of the four most
dominant eigenvalues with distinct and finite real part for different $N$ and
for $\beta\\!=\\!4$ as a representative of the SM phase. The data points
corresponding to system sizes larger than $N=350$ are not resulting from a
full diagonalization of the matrix but were obtained exploiting the sparseness
of the matrices (maximal 6 of the approximately $\\!N^{2}/2$ entries of every
row/column are nonzero), using a recursive algorithm, to obtain the dominant
part of the spectrum. The inset depicts the relaxation time scale $\tau_{r}$
and the lifetime of the metastable state $\tau_{l}$ as a function of $N$. In
b) the dashed, horizontal lines labeled as $\omega_{lc}$,$\omega_{ls}$
correspond to the LC frequency in the MF and to the imaginary part of the
linear stability matrix eigenvalue from Eq. (24), respectively. The mean
occupation density $\langle n_{i}\rangle$ as a function of time for both the
full (25b) and truncated (28) propagation [$i$=1 in c), d) and $i$=2 in e) and
f)] for the different network sizes $N=10^{2},10^{3}$. The dynamics for
$N\\!=\\!10^{3}$ was generated using the direct Gillespie method.

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 $\lambda_{1}$ smoothly converges
to the LC frequency $\omega_{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,
$\langle\bm{n}(t)\rangle\leavevmode\nobreak\ =\leavevmode\nobreak\
\sum_{\bm{N}}\bm{N}/N\,P_{\bm{N}}(t)$, using the full spectral decomposition
of the Markov generator in Eq. (25b) and the truncated one

$\displaystyle\bm{P}(t)\\!\overset{t\gg\tau_{m}}{\approx}\\!c_{0}\Phi^{R}_{0}+\mathrm{e}^{\lambda_{1}t}c_{1}\Phi^{R}_{1}+\mathrm{e}^{\lambda_{1}^{*}t}c_{1}^{*}\Phi^{R^{*}}_{1},$
(28)

for $N=10^{2},10^{3}$ at $\beta=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 $\beta=7$. In
contrast to Fig. 6a), here, $\operatorname{Re}(\lambda_{2})$ clearly converges
to a finite value with $\operatorname{Re}(\lambda_{1})$ quickly going to zero
already for small $N$. This is confirmed by the inset showing that $\tau_{l}$
and $\tau_{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).

Figure 7: Real a) and corresponding imaginary b) parts of the four most
dominant eigenvalues with distinct and finite real part for different $N$ and
for $\beta\\!=\\!7$ as a representative of the AM phase. The inset in a)
depicts the relaxation time scale $\tau_{r}$ and the lifetime of the
metastable state $\tau_{l}$ as a function of $N$. Furthermore, the mean
occupation density $\langle n_{i}\rangle$ as a function of time for both the
full (25b) and truncated (28) propagation [$i$=1 in c) and $i$=2 in d)] for
$N=10^{2}$ is depicted.

Thus, we conclude from the observations made in this section, that for
sufficiently large systems in the SM and AM phase at times
$\tau_{m}\\!\ll\\!t\\!\ll\\!\tau_{r}$, the relaxation dynamics is determined
by the metastable state associated with $\lambda_{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\gg\tau_{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\\!\sim\\!10^{3}$ via full
diagonalization of the propagator is computationally not feasible 222We
mention that the non-symmetric real matrix implies in general a complex
eigensystem which shall be determined with float precision. The amount of
random-access memory (RAM) required to diagonalize a matrix of dimension 80601
$\times$ 80601 corresponding to a system consisting of 400 units is about 312
GB. We restrict the diagonalization to that size and employ different
numerical methods for larger systems.. 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 Gillespie (1976, 1977), 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 $\langle\cdot\rangle$.

Figure 8 depicts the $\left\langle n_{2}\right\rangle\\!-\\!\left\langle
n_{1}\right\rangle$ plots generated with the Gillespie algorithm sampling over
$10^{6}$ trajectories for selected values of $\beta$ and for different system
sizes $N=10^{2},10^{4}$. Except for $\beta\\!=\\!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 ($\beta=4,5,6.1$) and AM phase ($\beta=7$). In the A phase
($\beta=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 $\beta=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, $\beta\approx\beta_{c_{2}}$,
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, $\omega_{lc}$, and the imaginary part of the crucial
eigenvalue, $\lambda_{1}$, as the second critical point, $\beta_{c_{2}}\approx
6.1068$ is approached [cf. Fig. 5c)].

Figure 8: Parametric plot of the mean occupation densities $\langle
n_{i}\rangle$ for different finite system sizes $N\\!=\\!10^{2}$ (blue solid
line), $N\\!=\\!10^{4}$ (green solid line), and the MF limit
$(N\\!=\\!\infty$, orange dotted line) at distinct values of $\beta$. In all
panels we initialize the system in the ground state with $n_{1}=1$ and sample
$10^{6}$ trajectories.

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),
$\bm{P}(0)\\!=\\!\sum_{i\neq 1}a_{i}\,\Phi^{R}_{i}$, excluding the mode
associated with the crucial eigenvalue pair $\lambda_{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 metastability would be removed from
the dynamics.

Figure 9: Joint probability distribution
$\mathfrak{P}\left\\{n_{1}(t),n_{2}(t)\right\\}$ at $\beta=4$ in a), b) and at
$\beta=7$ in c), d) for system sizes $N=10^{2},10^{3}$ at time
$t=\leavevmode\nobreak\ 20$. The plots were created using a grid of dimension
$101\times 101$ that specifies the set of initial conditions. For comparison,
in all plots the long-time MF solution (green solid line above, green closed
circles below) is overlaid.

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
$\mathfrak{P}\left\\{n_{1}(t=20),n_{2}(t=20)\right\\}$ for different system
sizes $N=10^{2},10^{4}$ and $\beta=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 $\beta=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
$\beta$ 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 $\beta$.
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
$\mathfrak{P}\left\\{n_{1}(t),n_{2}(t)\right\\}$ 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 333See supplementary media on
https://doi.org/10.6084/m9.figshare.5822097 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

$\displaystyle\left\langle e\right\rangle$
$\displaystyle=\sum\limits_{\alpha}E(\alpha)\,p_{\alpha}$ (29a)
$\displaystyle\left\langle s\right\rangle$
$\displaystyle=-\sum\limits_{\alpha}p_{\alpha}\ln p_{\alpha}\,.$ (29b)

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

$\displaystyle\mathrm{d}_{t}\langle
e\rangle=\sum\limits_{\alpha,\alpha^{\prime}}E(\alpha)w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}=\langle\dot{q}\rangle+\langle\dot{w}\rangle,$
(30)

naturally defines the microscopic first law of thermodynamics with the heat
and work rate given by

$\displaystyle\langle\dot{q}\rangle$
$\displaystyle=\sum\limits_{\alpha,\alpha^{\prime}}\left[E(\alpha)-f\,\Theta(\alpha,\alpha^{\prime})\right]w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}$
(31) $\displaystyle\langle\dot{w}\rangle$
$\displaystyle=\sum\limits_{\alpha,\alpha^{\prime}}f\;\Theta(\alpha,\alpha^{\prime})\;w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}\,,$
(32)

where the sign function $\Theta(\alpha,\alpha^{\prime})$ 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

$\displaystyle
q(\alpha,\alpha^{\prime})=-\frac{1}{\beta}\ln\frac{w_{\alpha\alpha^{\prime}}}{w_{\alpha^{\prime}\alpha}}.$
(33)

The system entropy change

$\displaystyle\mathrm{d}_{t}\langle
s\rangle=\langle\dot{s}_{e}\rangle+\langle\dot{\sigma}\rangle$ (34)

can be decomposed into the entropy flow from the bath to the system

$\displaystyle\langle\dot{s}_{e}\rangle=-\sum\limits_{\alpha,\alpha^{\prime}}w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}\ln\frac{w_{\alpha\alpha^{\prime}}}{w_{\alpha^{\prime}\alpha}}=\beta\langle\dot{q}\rangle,$
(35)

and the non-negative entropy production (EP) rate

$\displaystyle\langle\dot{\sigma}\rangle=\sum\limits_{\alpha,\alpha^{\prime}}w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}\,\ln\frac{w_{\alpha\alpha^{\prime}}p_{\alpha^{\prime}}}{w_{\alpha^{\prime}\alpha}p_{\alpha}}\geq
0.$ (36)

Equation (36) is the second law of thermodynamics and the inequality follows
straightforwardly from $\ln x\leq x-1$. The marginalization of the microscopic
probability $p_{\alpha}$ 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 Esposito (2012).
Defining $E_{\bm{N}}$ to be the internal energy of the system in the
macrostate $\bm{N}$, and applying the coarse-graining from Eq. (10) on the
expression for the internal energy in Eq. (29a), we obtain

$\displaystyle\mathrm{d}_{t}\langle e\rangle$
$\displaystyle=\sum\limits_{\bm{N},\bm{N}^{\prime}}\,\sum\limits_{\alpha^{\prime}\in\bm{N}^{\prime}}\,\sum\limits_{\alpha\in\bm{N}}E(\alpha)\,w_{\alpha\alpha^{\prime}}\,p_{\alpha^{\prime}}$
(37a)
$\displaystyle=\sum\limits_{\bm{N},\bm{N}^{\prime}}\tilde{W}_{\bm{N}\bm{N}^{\prime}}\,E(\bm{N})\\!\sum\limits_{\alpha^{\prime}\in\bm{N}^{\prime}}p_{\alpha^{\prime}}\sum\limits_{\alpha\in\bm{N}}1$
(37b)
$\displaystyle=\sum\limits_{\bm{N},\bm{N}^{\prime}}E(\bm{N})\,W_{\bm{N}\bm{N}^{\prime}}\,P_{\bm{N}^{\prime}}\equiv\mathrm{d}_{t}\langle
E\rangle.$ (37c)

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 coarse-grained as

$\displaystyle\langle\dot{Q}\rangle$
$\displaystyle=\sum\limits_{\bm{N},\bm{N}^{\prime}}\underbrace{\big{(}E(\bm{N})-f\,\Theta(\bm{N},\bm{N}^{\prime})\big{)}W_{\bm{N}\bm{N}^{\prime}}P_{\bm{N}^{\prime}}}_{Q(\bm{N},\bm{N}^{\prime})}$
(38a) $\displaystyle\langle\dot{W}\rangle$
$\displaystyle=\sum\limits_{\bm{N},\bm{N}^{\prime}}f\,\Theta(\bm{N},\bm{N}^{\prime})W_{\bm{N}\bm{N}^{\prime}}\,P_{\bm{N}^{\prime}}.$
(38b)

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

$\displaystyle\\!\\!\\!Q(\bm{N},\bm{N}^{\prime})$
$\displaystyle\\!=\\!-\frac{1}{\beta}\ln\frac{W_{\bm{N}\bm{N}^{\prime}}}{W_{\bm{N}^{\prime}\bm{N}}}\\!+\\!\frac{1}{\beta}\,\Delta
S^{int}(\bm{N},\bm{N}^{\prime})$ (39a)
$\displaystyle=-\frac{1}{\beta}\ln\frac{\tilde{W}_{\bm{N}\bm{N}^{\prime}}}{\tilde{W}_{\bm{N}^{\prime}\bm{N}}},$
(39b)

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)
Esposito (2012). We define the system entropy in the mesospace

$\displaystyle\left\langle S\right\rangle$
$\displaystyle=\sum\limits_{\bm{N}}P_{\bm{N}}\left(\Omega(\bm{N})-\ln
P_{\bm{N}}\right)\,,$ (40)

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

$\displaystyle\langle\dot{S}_{e}\rangle=-\sum\limits_{\bm{N},\bm{N}^{\prime}}W_{\bm{N}\bm{N}^{\prime}}\,P_{\bm{N}^{\prime}}\ln\frac{\tilde{W}_{\bm{N}\bm{N}^{\prime}}}{\tilde{W}_{\bm{N}^{\prime}\bm{N}}}=\beta\langle\dot{Q}\rangle,$
(41)

and the EP rate

$\displaystyle\langle\dot{\Sigma}\rangle=\sum\limits_{\bm{N},\bm{N}^{\prime}}W_{\bm{N}\bm{N}^{\prime}}\,P_{\bm{N}^{\prime}}\,\ln\frac{W_{\bm{N}\bm{N}^{\prime}}P_{\bm{N}^{\prime}}}{W_{\bm{N}^{\prime}\bm{N}}P_{\bm{N}}}\geq
0.$ (42)

The definitions in Eqs. (40),(42) are in general not coinciding with those
made at the microscopic level, _i.e._ $\langle S\rangle\neq\langle
s\rangle,\langle\Sigma\rangle\neq\langle\sigma\rangle$. The nonlinearity of
the system entropy and the EP [Eqs. (29b), (36)] in the microstate probability
$p_{\alpha}$ 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 Esposito (2012). For the special case of a stationary
probability distribution, $P^{s}_{\bm{N}}$, one can show (cf. appendix B) via
the spanning tree formula Schnakenberg (1976) that the microstates belonging
to the same macrostate are equally probable,
$p_{\alpha}=P_{\bm{N}}/\Omega(\bm{N})$. In the stationary limit, the entropies
in mesoscopic representation are therefore identical to those in microscopic
representation, _i.e._ $\langle S^{s}\rangle=\langle
s^{s}\rangle,\langle\Sigma^{s}\rangle=\langle\sigma^{s}\rangle$. For this
particular case, the second law

$\displaystyle\\!\\!\\!\\!\langle\dot{\Sigma^{s}}\rangle\\!=\\!\sum\limits_{\bm{N},\bm{N}^{\prime}}\\!\\!W_{\bm{N}\bm{N}^{\prime}}\,P^{s}_{\bm{N}^{\prime}}\,\ln\frac{\tilde{W}_{\bm{N}\bm{N}^{\prime}}}{\tilde{W}_{\bm{N}^{\prime}\bm{N}}}=-\langle\dot{S}_{e}^{s}\rangle\geq
0,$ (43)

boils down to the steady entropy flow $\langle\dot{S}_{e}^{s}\rangle$ being
equal to the magnitude of the steady EP rate $\langle\dot{\Sigma}^{s}\rangle$.
Using the non-positivity of the average stationary heat,
$\beta\langle\dot{Q}^{s}\rangle\leq 0$, we easily verify that
$\langle\dot{\Sigma}^{s}\rangle\geq 0$.

We now turn to the MF case and consistently define the first law in this limit

$\displaystyle\mathrm{d}_{t}\,\mathcal{E}=\sum\limits_{i,j}E_{i}\,k_{ij}\,\overline{n}_{j}=\dot{\mathcal{Q}}+\dot{\mathcal{W}},$
(44)

with the heat and work flux

$\displaystyle\dot{\mathcal{Q}}$
$\displaystyle=\sum\limits_{i,j}\big{(}E_{i}-f\,\Theta(i,j)\big{)}\,k_{ij}\,\overline{n}_{j}$
(45a) $\displaystyle\dot{\mathcal{W}}$
$\displaystyle=\sum\limits_{i,j}f\,\Theta(i,j)\,k_{ij}\,\overline{n}_{j}\,,$
(45b)

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

$\displaystyle\mathcal{S}$
$\displaystyle=-\sum\limits_{i}\overline{n}_{i}\ln\overline{n}_{i},$ (46)

which we split into the MF entropy flow

$\displaystyle\dot{\mathcal{S}}_{e}=-\sum\limits_{i,j}k_{ij}\,\overline{n}_{j}\ln\frac{k_{ij}}{k_{ji}}=\beta\dot{\mathcal{Q}},$
(47)

and the non-negative MF EP rate

$\displaystyle\dot{\mathcal{S}_{i}}=\sum\limits_{i,j}k_{ij}\,\overline{n}_{j}\,\ln\frac{k_{ij}\overline{n}_{j}}{k_{ji}\overline{n}_{i}}\geq
0.$ (48)

As the MF represents the asymptotic limit of the mesospace, it holds that all
the mesoscopic averages of the intensive observables $\langle X\rangle/N$ that
are consistent with the coarse-graining in Eq. (10) converge to the
corresponding observables $\mathcal{X}$ in the MF limit,
$\lim\limits_{N\to\infty}\tfrac{\langle\dot{X}\rangle}{N}=\dot{\mathcal{X}},$
with $X\\!\\!\leavevmode\nobreak\ =\leavevmode\nobreak\
\\!\\!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, $\overline{n}^{s}$, which yields for the second law in the MF
limit

$\displaystyle\dot{\mathcal{S}}_{i}^{s}=\sum\limits_{i,j}k_{ij}\,\overline{n}^{s}_{j}\ln\frac{k_{ij}}{k_{ji}}=-\dot{\mathcal{S}}_{e}^{s}\geq
0$ (49)

The non-negativity of the MF EP follows from the non-positivity 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, $\langle W\rangle>0$, and
dissipates it into the bath in the form of heat, $\langle Q\rangle<0$, for all
temperatures and system sizes.

Figure 10a) depicts the difference between the stationary work current of a
single unit, $\overline{W}_{1}=2\,\Gamma f\sinh\\!\left(f\beta/2\right)$, and
the asymptotic work current per unit in a network of size $N$,
$\overline{W}_{N}\\!\equiv\\!\frac{\langle W\rangle}{Nt}$ as a function of
$\beta$ for different $N$. The derivation of the single-unit stationary work
current, $\overline{W}_{1}$, is deferred to appendix C and given by Eq. (82).
The asymptotic work current, $\overline{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 $\beta_{c_{1,2}}$. The vicinities of these critical points
are magnified in the two insets. The phase transitions in the dissipated MF
work at $\beta_{c_{1}}$ and $\beta_{c_{2}}$ 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 finite-time
signatures of these nonequilibrium phase transitions at the bifurcation points
which blur out with decreasing system size.

In the high-temperature limit, $\beta\to 0$, the difference
$\Delta\overline{W}_{1N}\leavevmode\nobreak\ \equiv\leavevmode\nobreak\
\overline{W}_{1}\leavevmode\nobreak\ -\leavevmode\nobreak\ \overline{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\to 0$ as $N\to\infty$.). While for the MF this holds true
in the entire SA phase, for finite systems the range of $\beta$ 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, $\Delta\overline{W}_{1N}\\!>\\!0$. This work dissipation
gap, $\Delta\overline{W}_{1N}$, is a monotonically increasing function of
$\beta$ and becomes infinitely large in the low-temperature limit, since
$\Delta\overline{W}_{1N}/\overline{W}_{1}\to 1$ as $\beta\to\infty$. 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

$\displaystyle\lim_{\beta\to\infty}\overline{W}_{N}=\lim_{\beta\to\infty}\Gamma
f\left(\mathrm{e}^{\beta f}-1\right)\mathrm{e}^{-\frac{\beta(fN-Nu+u)}{2N}},$
(50)

which is subdominant to $\overline{W}_{1}$ [cf. Eq. (82)]

$\displaystyle\lim\limits_{\beta\to\infty}\frac{\Delta\overline{W}}{\overline{W}_{1}}=\lim\limits_{\beta\to\infty}1-\mathrm{e}^{\frac{\beta
u(N-1)}{2N}}=1.$ (51)

Figure 10: Panel a): Difference of the dissipated work for a single-unit,
$\overline{W}_{1}$, and for a unit in a network of size $N$,
$\overline{W}_{N}$, for inverse temperatures $\beta=0\ldots 10$. The time
$t=500$ is chosen to ensure that $\overline{W}_{N}$ has converged to its
(metastable) asymptotic value. Panel b): Difference of the dissipated work per
unit for networks of different size with $N<10^{4}$, for $\beta$ ranging from
0 to 10 and thus covering all three phases: Symmetric asynchronous phase (SA),
synchronous phase (S), non-symmetric asynchronous phase (NA) in the MF and
asynchronous phase (A), synchronous metastable phase (SM), and the
asynchronous metastable phase (AM) for finite metastable systems. As in panel
a), the time is $t=500$. The purple closed circles in panels a) and b)
represent the analytic expression given by Eqs. (51) and (52), respectively.
Panels c)–e): Plot of $\overline{W}_{N}$ for selected values at $\beta=2,4,7$
and system sizes $N\\!=\\!10^{2}$ (blue solid line) and $N\\!=\\!10^{4}$
(green solid line). This is the same data as the one underlying the blue solid
curve in plot b) but, for better visualization, the time $t$ is restricted
from 20 to 80. For comparison, the MF limit (orange dashed lines) is overlaid
in c)–e). In each plot all finite systems were simulated sampling $10^{6}$
trajectories.

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 decreased 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 system 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
$\Delta\overline{W}_{N,10^{4}}\equiv\overline{W}_{N}-\overline{W}_{10^{4}}<0$
while the opposite holds true in the AM phase, where
$\Delta\overline{W}_{N,10^{4}}>0$.

Again, we observe at the critical points significant changes in
$\Delta\overline{W}_{N,10^{4}}$: At the first critical point
$\Delta\overline{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^{\prime}=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),
$\Delta\overline{W}_{N,10^{4}}$ goes to zero in the high-temperature regime,
while in the low-temperature limit one obtains

$\displaystyle\lim\limits_{\beta\to\infty}\frac{\Delta\overline{W}_{N,10^{4}}}{\overline{W}_{\\!\\!N}}=\lim\limits_{\beta\to\infty}1-\mathrm{e}^{-\frac{\beta
u}{2}\left(\frac{1}{10^{4}}-\frac{1}{N}\right)}=1,$ (52)

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 $\beta$ 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 $\beta=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, $\Delta\overline{W}_{1\infty}$, as a
function of $\beta$ for different $f$. Again, $\Delta\overline{W}_{1\infty}=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
$\beta\geq\beta_{c_{1}}$ the second critical point is gradually shifting to
smaller $\beta$ [cf. Fig. 2] while the difference
$\Delta\overline{W}_{1\infty}$ 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.

Figure 11: Comparison between stationary single-unit work current
$\overline{W}_{1}$ and asymptotic MF work current
$\overline{W}_{\\!\\!\infty}$ as a function of $\beta$ for
$f=0.25,\,0.5,\,0.75,\,1.0$. The time $t=500$ is chosen such that the time-
averaged MF work has converged to its asymptotic value.

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 Cuetara and
Esposito (2015) 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}$ Vale and
Milligan (2000); Imparato (2015). 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 $\mathcal{W}_{1}^{s}$ and $\mathcal{W}_{2}^{s}$.
Substituting Eq. (23) into Eq. (49), yields the following decomposition of the
stationary EP in the MF limit

$\displaystyle\mathcal{S}_{i}^{s}=\mathcal{S}_{i}^{\mathcal{W}_{1}^{s}}+\mathcal{S}_{i}^{\mathcal{W}_{2}^{s}}\,,$
(53)

where
$\mathcal{S}_{i}^{\mathcal{W}_{k}^{s}}\\!=\\!\beta\,\mathcal{W}_{k}^{s},\,k\\!=\\!1,2$.
Based on Eq. (53), we use as an unambigious definition for the efficiency of
this work-to-work conversion (cf. Refs. Esposito  _et al._ (2012); Verley  _et
al._ (2014))

$\displaystyle\eta=-\frac{\mathcal{S}_{i}^{\mathcal{W}_{2}^{s}}}{\mathcal{S}_{i}^{\mathcal{W}_{1}^{s}}}=1-\frac{f}{f_{1}}.$
(54)

At equilibrium ($f\\!=\\!0$), the reversible limit, $\eta_{c}\\!=\\!1$ is
attained while out of equilibrium ($f\neq 0$) the efficiency is bounded,
$0\\!<\\!\eta\\!<\\!1$. Of particular interest is the efficiency at maximum
power (EMP) Curzon and Ahlborn (1975), which results from the optimization of
the stationary output power
$\mathcal{P}\equiv\partial\mathcal{W}_{2}^{s}/\partial t$ with respect to the
output force

$\displaystyle\eta^{*}=1-\left.\frac{f^{*}}{f_{1}}\right|_{f^{*}=\\!f_{1}\\!-\\!f_{2}^{*}},$
(55)

The maximization parameter $f^{*}$ is determined by the condition
$\partial\mathcal{P}/\partial f_{2}\\!=\\!0$, while fixing $f_{1}\\!=\\!1$ and
thus varying the total dissipation.

In the SA phase, $\beta<\beta_{c}$, the stationary power putput coincides with
the average work current of a single unit given by Eq. (82). 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 $\mathcal{P}$ as a function of $\beta$ and $f_{2}$ in a density plot.

Figure 12: Depiction of the output power a) as a function of the output force
$f_{2}$ and the inverse temperature $\beta$. 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 $\mathcal{P}^{*}$ is optimized with respect to
$f_{2}$ and in panel c) the associated EMP $\eta^{*}(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
$\beta<\beta_{c}$ [green lines] is overlaid with the numerical data in the
lower panels.

The white dashed lines indicate the critical points $\beta_{c_{1,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
$\beta$ 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

$\displaystyle\frac{\partial\mathcal{P}}{\partial f_{2}}$
$\displaystyle=\left[2+\beta f_{2}+\mathrm{e}^{\beta(1-f_{2})}(\beta
f_{2}-2)\right]=0,$ (56)

results in a transcendental equation that must be treated numerically. In the
low-temperature limit, the extremum condition

$\displaystyle\frac{\partial\mathcal{P}}{\partial f_{2}}$
$\displaystyle=\\!\mathrm{e}^{\frac{\beta}{2}(f_{2}-1)}\\!\left[\mathrm{e}^{\beta(1-f_{2})}(\beta
f_{2}\\!-\\!2)\\!+\\!(\beta f_{2}\\!+\\!2)\right]\\!=0,$ (57)

can not be satisfied for any $f_{2}$ compatible with the constraint
$\beta=\infty$.

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 $\beta$ approaches zero, the EMP takes the universal
linear-response value for tightly-coupled (only one net-current) systems,
$\eta^{*}=0.5\,\eta_{c}$ Van den Broeck (2005); Esposito  _et al._ (2009b).
This can be seen by expanding the expression for the stationary work current
in the SA phase given by Eq. (82) up to first order in $\beta$ which yields
the linear-response relation $J^{s}\\!\approx\\!L\,f$ with the Onsager
coefficient $L=\Gamma\,\beta$. Therefore, small products $\beta f$ correspond
to linear response in our model and lead to EMP values very close to 1/2. With
increasing $\beta$, 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-from-equilibrium 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 non-equilibrium phase
transitions for the performance of an assembly of nano-machines and suggests
synchronization as an operating mode faciliating very efficient energy-
conversion 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\to\infty$) 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\to\infty$).

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 Wood  _et al._ (2006a, b), new interesting spatiotemporal regimes
may emerge Escaff  _et al._ (2014).

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.

## Acknowledgments

T.H. thanks Hadrien Vroylandt and Artur Wachtel for insightful discussions. We
also mention that the simulations were carried out using the HPC facilities of
the University of Luxembourg Varrette  _et al._ (2014). This research was
supported by the National Research Fund, Luxembourg, in the frame of the AFR
PhD Grant 2016, No. 11271777 and by the European Research Council project
NanoThermo (ERC-2015-CoG Agreement No. 681456).

##  Appendix A Characterization of Hopf bifurcation

We shall in the following prove that the Hopf bifurcation observed in Sec. III
is supercritical, _i.e._ results in stable LCs. To characterize the LC close
to the bifurcation point, we consider the normal form of the Hopf bifurcation.
The procedure is detailed in Kuznetsov (1998).

At first, we transform the two-dimensional system in Eq. (21) into a single
equation

$\displaystyle\dot{z}$
$\displaystyle=\lambda(\Delta\beta)z+g(z,z^{*},\Delta\beta),$ (58)

where $z$ is a complex variable, $z^{*}$ its complex-conjugate,
$\Delta\beta=\beta-\beta_{c_{1}}$ gives the distance of the inverse
temperature to the critical inverse temperature of the Hopf bifurcation and
$g=O\big{(}\begin{Vmatrix}z\end{Vmatrix}\\!\textsuperscript{2}\big{)}$ is a
smooth function of $(z,z^{*},\Delta\beta)$.

Such a transformation is achieved by first finding the complex eigenvectors
$\bm{r}$ and $\bm{v}$ determined by

$\displaystyle\bm{J}(0)\bm{r}=\lambda(0)\bm{r},\quad\bm{J}(0)^{\top}\bm{v}=\lambda(0)^{*}\bm{v},$
(59)

where the real and non-symmetric Jacobian $\bm{J}$ resulting from the
linearization of Eq. (21) is evaluated at the bifurcation point
$\beta=\beta_{c_{1}}$, yielding

$\displaystyle\bm{r}$
$\displaystyle=\left(\frac{1}{2}(-1+\sqrt{3}\,\mathrm{i}),1\right)^{\top}$
(60) $\displaystyle\bm{v}$
$\displaystyle=\frac{1}{3-\sqrt{3}\,\mathrm{i}}\left(1+\sqrt{3}\,\mathrm{i},2\right)^{\top}.$
(61)

If $|\Delta\beta|$ is sufficiently small, the two-dimensional system from Eq.
(21) can be written as

$\displaystyle\dot{\bm{\overline{n}}}=\bm{J}(\Delta\beta)\bm{\overline{n}}+\bm{F}(\bm{\overline{n}},\Delta\beta),$
(62)

where $\bm{F}(\bm{\overline{n}},\Delta\beta)$ is a smooth vector function
whose components have Taylor expansions in $\bm{\overline{n}}$ starting with
at least quadratic terms,
$F_{1,2}$=$O\big{(}\begin{Vmatrix}\bm{\overline{n}}\end{Vmatrix}\\!\textsuperscript{2}\big{)}$.
Using Eq. (58) and the properties
$\langle\bm{v},\bm{r}\rangle=1$,$\langle\bm{v},\bm{r}^{*}\rangle=0$, one can
show that

$\displaystyle\\!\\!g(z,z^{*}\\!\\!,\Delta\beta)\\!=\\!\langle\bm{v}(\Delta\beta),\bm{F}(z\,\bm{r}(\Delta\beta)\\!+\\!z^{*}\bm{r}^{*}(\Delta\beta),\Delta\beta)\rangle.$
(63)

The function $g$ can be formally written as a Taylor series in the two complex
variables $z$ and $z^{*}$,

$\displaystyle g(z,z^{*},\Delta\beta)=\sum\limits_{k+l\geq
2}\frac{1}{k!l!}\frac{\partial^{k+l}}{\partial z^{k}\partial
z^{*^{l}}}g_{kl}(\Delta\beta)\,z^{k}z^{*^{l}},$ (64) with $\displaystyle
g_{kl}(\Delta\beta)\\!\\!=\\!\\!\\!\left.\langle\bm{v}(\Delta\beta),\\!\bm{F}(\\!z\bm{r}(\Delta\beta)\\!+\\!z^{\\!*}\bm{r}^{\\!*}(\Delta\beta),\Delta\beta)\rangle\right|_{z=0}\\!.$
(65)

Moreover, if the function $F(\bm{\overline{n}},\Delta\beta)$ from Eq. (62) is
represented as

$\displaystyle
F(\bm{x},0)=\frac{1}{2}B(\bm{x},\bm{x})+\frac{1}{6}C(\bm{x},\bm{x},\bm{x})+O\left(\begin{Vmatrix}x\end{Vmatrix}^{4}\right),$
(66)

where $B(\bm{x},\bm{y})$ and $C(\bm{x},\bm{y},\bm{u})$ are symmetric
multilinear vector functions of $\bm{x},\bm{y},\bm{u}\in\mathbb{R}^{2}$, it
follows that

$\displaystyle g_{20}$ $\displaystyle=\langle\bm{v},B(\bm{r},\bm{r})\rangle=0$
(67a) $\displaystyle g_{11}$
$\displaystyle=\langle\bm{v},B(\bm{r},\bm{r}^{*})\rangle=0$ (67b)
$\displaystyle g_{21}$
$\displaystyle=\langle\bm{v},C(\bm{r},\bm{r},\bm{r})\rangle.$ (67c)

In coordinates, one has for these vector functions

$\displaystyle B_{i}(\bm{x},\bm{y})$
$\displaystyle=\sum\limits_{j,k=1}^{2}\left.\frac{\partial^{2}F_{i}(\bm{\xi},0)}{\partial\xi_{j}\partial\xi_{k}}\right|_{\xi=0}\\!\\!\\!x_{j}\,y_{k},\;\;i\\!=\\!1,2$
(68) $\displaystyle C_{i}(\bm{x},\bm{y},\bm{u})$
$\displaystyle=\\!\\!\sum\limits_{j,k,l=1}^{2}\left.\frac{\partial^{3}F_{i}(\bm{\xi},0)}{\partial\xi_{j}\partial\xi_{k}\partial\xi_{l}}\right|_{\xi=0}\\!\\!\\!\\!\\!x_{j}y_{k}u_{l},\;i\\!=\\!1,2.$
(69)

With these expressions at hand, we can determine the first Lyapunov
coefficient $L_{1}$ as

$\displaystyle
L_{1}=\frac{1}{2\,\omega_{lc}^{2}}\,\operatorname{Re}\left(\mathrm{i}\,g_{20}\,g_{11}+\omega_{lc}\,g_{21}\right),$
(70)

where the eigenvalue of the Jacobian is decomposed as
$\lambda(\Delta\beta)$=$\sigma(\Delta\beta)$+$\mathrm{i}\,\omega(\Delta\beta)$
and

$\displaystyle\omega_{lc}=\left.\lambda(\Delta\beta)\right|_{\beta=\beta_{c_{1}}}=\Gamma\sqrt{3}\,\sinh\left(-\frac{3f}{2u}\right)$
(71)

is the LC frequency, $\omega_{lc}\equiv\omega(0)$, evaluated at the
bifurcation point $\Delta\beta$=0. For Eq. (70) to hold, the two requirements
$\omega(0)>0$ and $\sigma^{\prime}(0)<0$ must be met. From Eq. (71) and

$\displaystyle\sigma^{\prime}(0)=u\,\Gamma\,\cosh\left(\frac{3f}{2u}\right)\,$
(72)

it follows that this is only true for attractive interactions, $u<0$.
Collecting results, we finally arrive at

$\displaystyle L_{1}=-\frac{81}{2}\Gamma\cosh\left(\frac{3f}{2u}\right),$ (73)

which is negative for any $u<0$, hence for attractive interactions stable LCs
emerge at the bifurcation point $\beta_{c_{1}}$ as asserted above.

##  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,
$\mathcal{T}^{\prime}(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
$\mathcal{T}^{\prime}(G)$ is connected and contains no circuits. We introduce
the notation $\mathcal{A}(\mathcal{T}^{\prime(\mu)}_{\alpha}(G))$ referring to
the $\mu$th spanning tree rooted in $\alpha$, that is a tree whose branches
are pointing towards the vertex $\alpha$. The spanning tree formula states
Schnakenberg (1976)

$\displaystyle
p^{s}_{\alpha}\\!\\!=\\!\frac{\sum\limits_{\mu}\mathcal{A}(\mathcal{T}^{\prime(\mu)}_{\alpha}(G))}{\sum\limits_{\alpha}\sum\limits_{\mu}\mathcal{A}(\mathcal{T}^{\prime(\mu)}_{\alpha}(G))}\\!=\\!\frac{\sum\limits_{\mathcal{T}^{\prime}_{\alpha}(G)}\prod\limits_{\begin{subarray}{c}\text{s.t.
current}\\\ \text{ is directed to
}\alpha\end{subarray}}w_{\alpha^{\prime}\alpha^{\prime\prime}}}{\sum\limits_{\alpha}\sum\limits_{\mathcal{T}^{\prime}_{\alpha}(G)}\prod\limits_{\begin{subarray}{c}\text{s.t.
current}\\\ \text{ is directed to
}\alpha\end{subarray}}w_{\alpha^{\prime}\alpha^{\prime\prime}}}\,.$ (74)

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 $\alpha$ is constant for all $\alpha$ inside the same macrostate.
Thus, at steady state, all microstates constituting the same macrostate

$\displaystyle p^{s}_{\alpha}\\!$
$\displaystyle=\\!\frac{\\!\\!\\!\\!\sum\limits_{\mathcal{T}^{\prime}_{\alpha\in\bm{N}}(G)}\prod\limits_{\begin{subarray}{c}\text{s.t.
current}\\\ \text{ is directed to
}\alpha\in\bm{N}\end{subarray}}\left.w_{\alpha\alpha^{\prime}}\right|_{\begin{subarray}{c}\alpha\in\bm{N}\\\
\alpha^{\prime}\in\bm{N}^{\prime}\end{subarray}}}{\sum\limits_{\alpha\in\bm{N}}\sum\limits_{\mathcal{T}^{\prime}_{\alpha(N)}(G)}\\!\\!\prod\limits_{\begin{subarray}{c}\text{s.t.
current}\\\ \text{ is directed to
}\alpha\in\bm{N}\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left.w_{\alpha\alpha^{\prime}}\right|_{\begin{subarray}{c}\alpha\in\bm{N}\\\
\alpha^{\prime}\in\bm{N}^{\prime}\end{subarray}}}\\!=\\!\mathrm{const},$ (75)

are equally probable and hence

$\displaystyle p^{s}_{\alpha}$
$\displaystyle=\frac{P^{s}_{\bm{N}}}{\Omega(\bm{N})}\,,$ (76)

where $\Omega(\bm{N})$ is the number of microstates forming the macrostate
$\bm{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

$\displaystyle P_{i}=\sum\limits_{i,j}W_{ij}\,P_{j}\,,$ (77)

where $P$ is the (macro-)probability to find the unit in the single state $i$
with the transition rates

$\displaystyle
W_{ij}=\mathrm{e}^{-\frac{\beta}{2}\left(\epsilon_{i}-\epsilon_{j}+\Theta(i,j)\,f\right)}\,,$
(78)

with the sign function $\Theta(i,j)$ as defined in Eq. (23) ensuring the
validity of local detailed balance. The steady-state work current reads

$\displaystyle\langle\dot{W}^{s}\rangle$
$\displaystyle=f\,\sum\limits_{i,j}\Theta(i,j)\,W_{ij}\,P^{s}_{\j}\,.$ (79)

Using the spanning tree formula from Eq. (74), one obtains for the stationary
probabilities

$\displaystyle P_{1}^{s}$ $\displaystyle=\frac{a_{1}}{a_{1}+a_{2}+a_{3}},\quad
P_{2}^{s}=\frac{a_{2}}{a_{1}+a_{2}+a_{3}}\,,$ (80)

where

$\displaystyle a_{1}$
$\displaystyle=W_{13}W_{12}\\!+\\!W_{12}W_{23}\\!+\\!W_{13}W_{32}$ (81a)
$\displaystyle a_{2}$
$\displaystyle=W_{23}W_{31}\\!+\\!W_{21}W_{13}\\!+\\!W_{23}W_{21}$ (81b)
$\displaystyle a_{3}$
$\displaystyle=W_{31}W_{12}\\!+\\!W_{32}W_{21}\\!+\\!W_{32}W_{31},$ (81c)

For a flat energy landscape, $\epsilon_{i}$=$\operatorname{const}$, we indeed
find that the symmetric stationary solution $P_{i}=1/3$ is independent of
$\beta$ and $f$ like in the MF limit. Next, the stationary work current is
given by

$\displaystyle\langle\dot{W}^{s}\rangle$
$\displaystyle=3f\,\frac{W_{13}W_{21}W_{32}-W_{31}W_{12}W_{23}}{W_{12}(W_{13}+W_{23}+W_{31})+W_{13}(W_{21}+W_{32})+(W_{21}+W_{31})(W_{13}+W_{23}+W_{31})}$
(82)

that simplifies to $\langle\dot{\mathcal{W}}^{s}\rangle$=$2\,\Gamma
f\sinh\left(f\beta/2\right)$ (see Sec. VII).

## References

  * Strogatz (2001) S. H. Strogatz,  _Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering_ (Westview Press, 2001).
  * Guckenheimer and Holmes (1983) John Guckenheimer and P. J. Holmes,  _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_ (Springer-Verlag New York, 1983).
  * Seifert (2012) Udo Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Reports on Progress in Physics 75, 126001 (2012).
  * den Broeck and Esposito (2015) C. Van den Broeck and M. Esposito, “Ensemble and trajectory thermodynamics: A brief introduction,” Physica A: Statistical Mechanics and its Applications 418, 6 – 16 (2015), proceedings of the 13th International Summer School on Fundamental Problems in Statistical Physics.
  * Zhang  _et al._ (2012) Xue-Juan Zhang, Hong Qian, and Min Qian, “Stochastic theory of nonequilibrium steady states and its applications. part i,” Physics Reports 510, 1 – 86 (2012).
  * Rao and Esposito (2018) Riccardo Rao and Massimiliano Esposito, “Conservation laws shape dissipation,” New Journal of Physics 20, 023007 (2018).
  * Pekola (2015) Jukka P Pekola, “Towards quantum thermodynamics in electronic circuits,” Nature Physics 11, 118–123 (2015).
  * Ciliberto (2017) S. Ciliberto, “Experiments in stochastic thermodynamics: Short history and perspectives,” Phys. Rev. X 7, 021051 (2017).
  * Esposito  _et al._ (2009a) M. Esposito, K. Lindenberg, and C. Van den Broeck, “Thermoelectric efficiency at maximum power in a quantum dot,” EPL (Europhysics Letters) 85, 60010 (2009a).
  * Cuetara and Esposito (2015) Gregory Bulnes Cuetara and Massimiliano Esposito, “Double quantum dot coupled to a quantum point contact: a stochastic thermodynamics approach,” New Journal of Physics 17, 095005 (2015).
  * Rutten  _et al._ (2009) B. Rutten, M. Esposito, and B. Cleuren, “Reaching optimal efficiencies using nanosized photoelectric devices,” Phys. Rev. B 80, 235122 (2009).
  * Ge  _et al._ (2012) Hao Ge, Min Qian, and Hong Qian, “Stochastic theory of nonequilibrium steady states. part ii: Applications in chemical biophysics,” Physics Reports 510, 87 – 118 (2012).
  * Gaspard and Gerritsma (2007) P. Gaspard and E. Gerritsma, “The stochastic chemomechanics of the f1-atpase molecular motor,” Journal of Theoretical Biology 247, 672 – 686 (2007).
  * Lau  _et al._ (2007) A. W. C. Lau, D. Lacoste, and K. Mallick, “Nonequilibrium fluctuations and mechanochemical couplings of a molecular motor,” Phys. Rev. Lett. 99, 158102 (2007).
  * Seifert (2011) U. Seifert, “Stochastic thermodynamics of single enzymes and molecular motors,” The European Physical Journal E 34, 26 (2011).
  * Altaner  _et al._ (2015) Bernhard Altaner, Artur Wachtel, and Jürgen Vollmer, “Fluctuating currents in stochastic thermodynamics. ii. energy conversion and nonequilibrium response in kinesin models,” Phys. Rev. E 92, 042133 (2015).
  * Shiraishi  _et al._ (2016) Naoto Shiraishi, Keiji Saito, and Hal Tasaki, “Universal trade-off relation between power and efficiency for heat engines,” Phys. Rev. Lett. 117, 190601 (2016).
  * Proesmans  _et al._ (2016) Karel Proesmans, Bart Cleuren, and Christian Van den Broeck, “Power-efficiency-dissipation relations in linear thermodynamics,” Phys. Rev. Lett. 116, 220601 (2016).
  * Pietzonka and Seifert (2018) Patrick Pietzonka and Udo Seifert, “Universal trade-off between power, efficiency, and constancy in steady-state heat engines,” Phys. Rev. Lett. 120, 190602 (2018).
  * Lee and Park (2017) Jae Sung Lee and Hyunggyu Park, “Carnot efficiency is reachable in an irreversible process,” Scientific Reports 7, 2045–2322 (2017).
  * Polettini and Esposito (2017) Matteo Polettini and Massimiliano Esposito, “Carnot efficiency at divergent power output,” EPL (Europhysics Letters) 118, 40003 (2017).
  * Vroylandt  _et al._ (2017) Hadrien Vroylandt, Massimiliano Esposito, and Gatien Verley, “Collective effects enhancing power and efficiency,” EPL (Europhysics Letters) 120, 30009 (2017).
  * Shiraishi (2017) Naoto Shiraishi, “Stationary engines in and beyond the linear response regime at the carnot efficiency,” Phys. Rev. E 95, 052128 (2017).
  * Golubeva and Imparato (2012) N. Golubeva and A. Imparato, “Efficiency at maximum power of interacting molecular machines,” Phys. Rev. Lett. 109, 190602 (2012).
  * Golubeva and Imparato (2013) N. Golubeva and A. Imparato, “Maximum power operation of interacting molecular motors,” Phys. Rev. E 88, 012114 (2013).
  * Imparato (2015) Alberto Imparato, “Stochastic thermodynamics in many-particle systems,” New Journal of Physics 17, 125004 (2015).
  * Sasa (2015) S. Sasa, “Collective dynamics from stochastic thermodynamics,” New Journal of Physics 17, 045024 (2015).
  * Pikovsky  _et al._ (2003) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_, Vol. 12 (Cambridge University Press, 2003).
  * Huygens (1895) Letters to de Sluse (letters; no. 1333 of 24 February 1665 no. 1335 of 26 February 1665 no. 1345 of 6 March 1665) Huygens, C., (Societe Hollandaise Des Sciences, Martinus Nijho, La Haye, 1895).
  * Strogatz (2004) S. Strogatz,  _Sync: The emerging science of spontaneous order_ , 1st ed. (Penguin UK, 2004).
  * Winfree (1967) Arthur T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” Journal of Theoretical Biology 16, 15 – 42 (1967).
  * Kuramoto (2003) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_, Dover Books on Chemistry Series (Dover Publications, 2003).
  * Strogatz (2000) Steven H. Strogatz, “From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators,” Physica D: Nonlinear Phenomena 143, 1 – 20 (2000).
  * Acebrón  _et al._ (2005) Juan A. Acebrón, L. L. Bonilla, Conrad J. Pérez Vicente, Félix Ritort, and Renato Spigler, “The kuramoto model: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys. 77, 137–185 (2005).
  * Gupta  _et al._ (2014) Shamik Gupta, Alessandro Campa, and Stefano Ruffo, “Kuramoto model of synchronization: equilibrium and nonequilibrium aspects,” Journal of Statistical Mechanics: Theory and Experiment 2014, R08001 (2014).
  * Sakaguchi  _et al._ (1987) Hidetsugu Sakaguchi, Shigeru Shinomoto, and Yoshiki Kuramoto, “Local and grobal self-entrainments in oscillator lattices,” Progress of Theoretical Physics 77, 1005–1010 (1987).
  * Daido (1988) Hiroaki Daido, “Lower critical dimension for populations of oscillators with randomly distributed frequencies: A renormalization-group analysis,” Phys. Rev. Lett. 61, 231–234 (1988).
  * Walgraef  _et al._ (1983) D. Walgraef, G. Dewel, and P. Borckmans, “Chemical waves in a two‐dimensional oscillating system,” The Journal of Chemical Physics 78, 3043–3051 (1983).
  * Izumida  _et al._ (2016) Yuki Izumida, Hiroshi Kori, and Udo Seifert, “Energetics of synchronization in coupled oscillators rotating on circular trajectories,” Phys. Rev. E 94, 052221 (2016).
  * Risler  _et al._ (2004) Thomas Risler, Jacques Prost, and Frank Jülicher, “Universal critical behavior of noisy coupled oscillators,” Phys. Rev. Lett. 93, 175702 (2004).
  * Risler  _et al._ (2005) Thomas Risler, Jacques Prost, and Frank Jülicher, “Universal critical behavior of noisy coupled oscillators: A renormalization group study,” Phys. Rev. E 72, 016130 (2005).
  * Wood  _et al._ (2006a) Kevin Wood, C. Van den Broeck, R. Kawai, and Katja Lindenberg, “Universality of synchrony: Critical behavior in a discrete model of stochastic phase-coupled oscillators,” Phys. Rev. Lett. 96, 145701 (2006a).
  * Wood  _et al._ (2006b) Kevin Wood, C. Van den Broeck, R. Kawai, and Katja Lindenberg, “Critical behavior and synchronization of discrete stochastic phase-coupled oscillators,” Phys. Rev. E 74, 031113 (2006b).
  * Wood  _et al._ (2007) Kevin Wood, C. Van den Broeck, R. Kawai, and Katja Lindenberg, “Continuous and discontinuous phase transitions and partial synchronization in stochastic three-state oscillators,” Phys. Rev. E 76, 041132 (2007).
  * Van Kampen (2007) N.G. Van Kampen,  _Stochastic Processes in Physics and Chemistry_, 3rd ed. (North Holland, 2007).
  * Gaspard (1996) Pierre Gaspard, “Hydrodynamic modes as singular eigenstates of the liouvillian dynamics: Deterministic diffusion,” Phys. Rev. E 53, 4379–4401 (1996).
  * Gaspard (1998) P. Gaspard,  _Chaos, Scattering, and Statistical Mechanics_ (Cambridge University Press, 1998).
  * Gaveau and Schulman (1998) Bernard Gaveau and L. S. Schulman, “Theory of nonequilibrium first-order phase transitions for stochastic dynamics,” Journal of Mathematical Physics 39, 1517–1533 (1998).
  * Biroli and Kurchan (2001) Giulio Biroli and Jorge Kurchan, “Metastable states in glassy systems,” Phys. Rev. E 64, 016101 (2001).
  * Esposito and Gaspard (2005a) M. Esposito and P. Gaspard, “Emergence of diffusion in finite quantum systems,” Phys. Rev. B 71, 214302 (2005a).
  * Esposito and Gaspard (2005b) Massimiliano Esposito and Pierre Gaspard, “Exactly Solvable Model of Quantum Diffusion,” J. Stat. Phys. 121, 463–496 (2005b).
  * Macieszczak  _et al._ (2016) Katarzyna Macieszczak, Madalin Guta, Igor Lesanovsky, and Juan P. Garrahan, “Towards a theory of metastability in open quantum dynamics,” Phys. Rev. Lett. 116, 240404 (2016).
  * Rose  _et al._ (2016) Dominic C. Rose, Katarzyna Macieszczak, Igor Lesanovsky, and Juan P. Garrahan, “Metastability in an open quantum ising model,” Phys. Rev. E 94, 052132 (2016).
  * Campàs  _et al._ (2006) O. Campàs, Y. Kafri, K. B. Zeldovich, J. Casademunt, and J.-F. Joanny, “Collective dynamics of interacting molecular motors,” Phys. Rev. Lett. 97, 038101 (2006).
  * Hill (1977) T. L. Hill, “Theoretical study of the effect of enzyme-enzyme interactions on steady-state enzyme kinetics,” PNAS 74, 3632–3636 (1977).
  * Hill and Chen (1981) T. L. Hill and Y. D. Chen, “Three-state, steady-state Ising systems: Monte Carlo and Bragg-Williams treatments,” PNAS 78, 4–8 (1981).
  * Raymo  _et al._ (2003) Françisco M. Raymo, Silvia Giordani, Andrew J. P. White, and David J. Williams, “Digital Processing with a Three-State Molecular Switch,” J. Org. Chem. 68, 4158–4169 (2003).
  * Simão  _et al._ (2011) Cláudia Simão, Marta Mas-Torrent, Javier Casado-Montenegro, Francisco Otón, Jaume Veciana, and Concepció Rovira, “A Three-State Surface-Confined Molecular Switch with Multiple Channel Outputs,” J. Am. Chem. Soc. 133, 13256–13259 (2011).
  * Note (1) Note that a state does not correspond to a folding of the two triangles but a rotation of one of the two planes about $180^{\circ}$ and subsequent overlapping of the two layers.
  * Keener (1981) Jamse P Keener, “Infinite period bifurcation and global bifurcation branches,” SIAM Journal on Applied Mathematics 41, 127–144 (1981).
  * Note (2) We mention that the non-symmetric real matrix implies in general a complex eigensystem which shall be determined with float precision. The amount of random-access memory (RAM) required to diagonalize a matrix of dimension 80601 $\times$ 80601 corresponding to a system consisting of 400 units is about 312 GB. We restrict the diagonalization to that size and employ different numerical methods for larger systems.
  * Gillespie (1976) Daniel T Gillespie, “A general method for numerically simulating the stochastic time evolution of coupled chemical reactions,” J. Comp. Phys. 22, 403–434 (1976).
  * Gillespie (1977) Daniel T Gillespie, “Exact stochastic simulation of coupled chemical reactions,” J. Phys. Chem. 81, 2340–2361 (1977).
  * Note (3) See supplementary media on https://doi.org/10.6084/m9.figshare.5822097.
  * Esposito (2012) Massimiliano Esposito, “Stochastic thermodynamics under coarse graining,” Phys. Rev. E 85, 041125 (2012).
  * Schnakenberg (1976) J. Schnakenberg, “Network theory of microscopic and macroscopic behavior of master equation systems,” Rev. Mod. Phys. 48, 571–585 (1976).
  * Vale and Milligan (2000) Ronald D. Vale and Ronald A. Milligan, “The way things move: Looking under the hood of molecular motor proteins,” Science 288, 88–95 (2000).
  * Esposito  _et al._ (2012) Massimiliano Esposito, Niraj Kumar, Katja Lindenberg, and Christian Van den Broeck, “Stochastically driven single-level quantum dot: A nanoscale finite-time thermodynamic machine and its various operational modes,” Phys. Rev. E 85, 031117 (2012).
  * Verley  _et al._ (2014) Gatien Verley, Tim Willaert, Christian Van den Broeck, and Massimiliano Esposito, “The unlikely carnot efficiency,” Nature Communications 5, 4721 (2014).
  * Curzon and Ahlborn (1975) F. L. Curzon and B. Ahlborn, “Efficiency of a Carnot engine at maximum power output,” American Journal of Physics 43, 22–24 (1975).
  * Van den Broeck (2005) C. Van den Broeck, “Thermodynamic efficiency at maximum power,” Phys. Rev. Lett. 95, 190602 (2005).
  * Esposito  _et al._ (2009b) Massimiliano Esposito, Katja Lindenberg, and Christian Van den Broeck, “Universality of efficiency at maximum power,” Phys. Rev. Lett. 102, 130602 (2009b).
  * Escaff  _et al._ (2014) Daniel Escaff, Italo’Ivo Lima Dias Pinto, and Katja Lindenberg, “Arrays of stochastic oscillators: Nonlocal coupling, clustering, and wave formation,” Phys. Rev. E 90, 052111 (2014).
  * Varrette  _et al._ (2014) S. Varrette, P. Bouvry, H. Cartiaux, and F. Georgatos, “Management of an academic hpc cluster: The ul experience,” in  _Proc. of the 2014 Intl. Conf. on High Performance Computing & Simulation (HPCS 2014)_ (IEEE, Bologna, Italy, 2014) pp. 959–967.
  * Kuznetsov (1998) Y. A. Kuznetsov,  _Elements of applied bifurcation theory_, 2nd ed. (Springer-Verlag, 1998).