Stochastic thermodynamics of all-to-all interacting many-body systems

We provide a stochastic thermodynamic description across scales for N identical units with all-to-all interactions that are driven away from equilibrium by different reservoirs and external forces. We start at the microscopic level with Poisson rates describing transitions between many-body states. We then identify an exact coarse graining leading to a mesoscopic description in terms of Poisson transitions between systems occupations. We also study macroscopic fluctuations using the Martin-Siggia-Rose formalism and large deviation theory. In the macroscopic limit ($N \to \infty$), we derive an exact nonlinear (mean-field) rate equation describing the deterministic dynamics of the most likely occupations. Thermodynamic consistency, in particular the detailed fluctuation theorem, is demonstrated across microscopic, mesoscopic and macroscopic scales. The emergent notion of entropy at different scales is also outlined. Macroscopic fluctuations are calculated semi-analytically in an out-of-equilibrium Ising model. Our work provides a powerful framework to study thermodynamics of nonequilibrium phase transitions.


I. INTRODUCTION
Interacting many body systems can give rise to a very rich variety of emergent behaviors such as phase transitions. At equilibrium, their thermodynamic properties have been the object of intensive studies and are nowadays well understood [1][2][3][4], see also [5] for a more philosophical perspective. When driven out-of-equilibrium, these systems are known to give rise to complex dynamical behaviors [6][7][8][9][10][11][12]. While most of the works are focused on their ensemble averaged description, in recent years progress was also made in characterizing their fluctuations [13][14][15][16][17]. However, little is known about their thermodynamic description. For instance, thermodynamics of nonequilibrium phase transitions started to be explored only recently [18][19][20][21][22][23][24][25][26][27][28][29][30]. There is a pressing need to develop methodologies to study thermodynamic quantities such as heat work and dissipation, not only at the average but also at the fluctuation level. To do so one has to start from stochastic thermodynamics that has proven instrumental to systematically infer the thermodynamics of small systems that can be driven arbitrarily far from equilibrium [31][32][33][34][35][36]. This theory consistently builds thermodynamics on top of a Markov dynamics (e.g. master equations [37] or Fokker-Planck equations [38]) describing open systems interacting with their surrounding. Its predictions have been experimentally validated in a broad range of fields ranging from electronics to single molecules and Brownian particles [39,40]. It has been particularly successful in studying the performance of small energy converters operating far-from-equilibrium and their power-efficiency trade-off [31,[41][42][43][44][45]. Until now, most of the focus has been on systems with finite phase space or few particle systems. However there are exceptions. Interacting systems have started to be considered in the context of energy conversion to asses whether they can trigger synergies in large ensembles of interacting energy converters. Beside few works such as [18,19], most other studies are restricted to mean-field treatments [20-24, 46, 47]. Another exception are chemical reaction networks which provide an interesting class of interacting systems. Indeed, while molecules in ideal solution are by definition noninteracting from an energetic standpoint, the stoichiometry of non-unimolecular reactions creates correlations amongst molecular species which generate entropic interactions. In the macroscopic limit, the mean field dynamics is exact and nonlinear [48][49][50] and can give rise to all sorts of complex behaviors [51]. The thermodynamics of chemical reaction networks has started to raise some attention in recent years [52][53][54][55][56].
The main achievement of this paper is to provide a consistent nonequilibrium thermodynamic description across scales of many body systems with all-to-all interactions. We do so by considering N identical units with all-to-all (or infinite range) interactions. Each unit is composed of q discrete states and undergoes transitions caused by one or more reservoirs. It may also be driven by an external force. The thermodynamics of this open many-body system is formulated at the ensemble averaged and fluctuating level, for finite N as well as in the macroscopic limit N → ∞. At the microscopic level, the system is characterized by microstates which correspond to the many-body states (i.e. they define the state of each of the units). Poisson rates describe the transitions between the microstates triggered by the reservoirs. These rates satisfy local detail balance, i.e. their log-ratio is the entropy change in the reservoir caused by the transition [57]. It implicitly assumes that the system is weakly coupled to reservoirs which instantaneously relax back to equilibrium after an exchange with the system. By linking the stochastic dynamics with the physics, this crucial property ensures a consistent nonequilibrium thermodynamics description of the system, in particular a detailed fluctuation theorem and an ensuing second law at the ensemble averaged level. The entropy of a state is given by minus the logarithm of the probability to find the system in that state and the ensemble averaged entropy is the corresponding Shannon entropy. Because we assume all units to be identical in they they interact with each other and in the way they interact with the reservoirs, we show that the microscopic stochastic dynamics can be exactly coarse grained to a mesoscopic level, where each system state specifies the unit occupations (i.e. the exact number of units which are in each of the unit states). The mesoscopic rates describing transition between occupations satisfy a local detailed balance. At this level, the entropy of a state is given by minus the logarithm of the probability to find the system in that state plus the internal entropy given by the logarithm of the number of microstates inside a mesostate, reflecting the fact that the units are energetically indistinguishable. We demonstrate that stochastic thermodynamics is invariant under this exact coarse-graining of the stochastic dynamics, if one considers initial conditions which are uniform within each mesostate, or for systems in stationary states. We then consider the macroscopic limit (N → ∞). Using a path integral representation of the stochastic dynamics (Martin-Siggia-Rose formalism), we also prove that the macroscopic fluctuations (i.e. the fluctuations that scale exponentially with the system size) satisfy a detailed fluctuation theorem and are thus thermodynamically consistent. We show via the path-integral representation that the stochastic dynamics exactly reduces to a mean-field rate equations with nonlinear rates and for deterministic variables corresponding to the most likely values of the occupation of each unit state. Remarkably, the nonlinear rates still satisfy local detailed balance and the entropy of each deterministic occupation is given by minus their logarithm. The entropy is thus a Shannon entropy for deterministic variables exclusively arising from the entropy inside the mesostates and not from the probability distribution to be on a mesostates. Indeed, this latter narrows down around its single or multiple (in case of phase transition) most likely values and gives rise to a vanishing stochastic entropy. We finally use our methodology to calculate macroscopic fluctuations in a semi-analytically solvable Ising model in contact with two reservoirs and displaying a nonequilibrium phase transition. The plan of the paper is as follows. First, in Sec. II, the many-body model is introduced and the stochastic dynamics is formulated. Moreover, the exact coarse-graining scheme is presented and the asymptotic mean-field equations are derived. Next, in Sec. III, using the formalism of stochastic thermodynamics and Martin-Siggia-Rose, the fluctuating thermodynamic quantities are formulated at different scales and the conditions under which they are preserved across these scales are identified. These theoretical results are illustrated via a semi-analytically solvable Ising model. We conclude with a summary and perspectives in Sec. VI.

A. Microscopic Description
We consider a system that consists of N all-to-all interacting identical and classical units that consist of q states i with energies i (λ t ) that are varying in time according to a known protocol λ t of an external driving. The system is coupled with multiple heat reservoirs ν = 1, 2, . . . , L at inverse temperatures β (ν) . Each unit is assumed to be fully connected, i.e. any state of a given unit can be reached within a finite number of steps from all other states of that unit, so that the global system is irreducible. Moreover, we suppose that all units are subjected to generic nonconservative forces f (ν) ij . Depending on whether a transition is aligned with or acting against the nonconservative force, the latter fosters or represses the transition from state j to i. For generality, the force is assumed to be different depending on which heat reservoir ν the system is exchanging energy with during the transition from j to i. Until explicitly states otherwise, we will take N to be finite in the following.
The many-body system is unambiguously characterized by a microstate α = (α 1 , . . . , α i , . . . , α N ), α i = 1, 2, . . . , q. (1) The system energy consists of the state occupation of the units and the interactions between them. For all-to-all interactions, we readily determine the energy of the system in a microstate α as follows, where u i (λ t )/N and u ij (λ t )/N denote the pair potential of units occupying the same or different single-unit states, respectively. These interactions can be tuned by an external driving according to a known protocol λ t , hence λ t = λ t , λ t . Moreover, N i (α) refers to the number of units N i occupying the single-unit state i for a given microstate α.
The stochastic jump process is governed by an irreducible Markovian master equation which describes the time evolution of the microscopic probability p α for the system to be in the microstate α as follows, with the microscopic rates w αα (λ t ) for transitions from α to α that in general depend on the current value of the driving parameter λ t . We note that probability conservation is ensured by the stochastic property of the transition rate matrix, α w αα (λ t ) = 0. The transition from α to α is induced by one of the L heat reservoirs, thus Here, for simplicity we assume that the transition rates are additive in the reservoirs ν. A more general treatment can be made following the procedure described in Ref. [57]. The microscopic transition rates that specify the heat reservoir satisfy the microscopic local detailed balance condition separately, which in turn ensures the thermodynamic consistency of the system. Here, f αα is the element of the nonconservative force vector f (ν) that is equal to f (ν) ij , if the microscopic transition from α → α corresponds to a single-unit transition from j → i. If the transition rates are kept constant, λ t = λ, the dynamics will relax into a unique stationary state, ∂ t p s α (λ) = 0. If furthermore all heat reservoirs have the same inverse temperature, β (ν) = β ∀ν, and the nonconservative forces vanish, f (ν) = 0 ∀ν, the stationary distribution coincides with the equilibrium one which satisfies the microscopic detailed balance condition, The local detailed balance (5) implies that the microscopic equilibrium distribution assumes the canonical form, with the microscopic equilibrium free energy

B. Mesoscopic Description
The microscopic state space grows exponentially with the number of units, ||α|| = q N . Yet, the complexity of the system can be significantly reduced. First, we note that due to the all-to-all interactions, there are equi-energetic microstates that are characterized by the same values for the occupation numbers N i . Next, we assume that the units are not only indistinguishable energetically [asymmetric part of the microscopic transition rates (5)] but also kinetically [symmetric part of the microscopic transition rates (4)] because they are all coupled in the same way to the reservoirs. As a result, the microscopic transition rates do not depend on the detailed pair of microstates that they connect but only on the pair of mesostate N ≡ (N 1 , N 2 , . . . , N q ) that they connect.
Consequentially, the microscopic dynamics can be marginalized into a mesoscopic one, where the mesostate N now identifies the state of the system. We denote by α N the equienergetic microstates α inside a mesostate N , that is microstates for which the relation holds. The number Ω N of microstates which belong to a mesostate is given by We introduce the mesoscopic probability to observe the mesostate N The conditional probability to find the system in a microstate α N that belongs to that mesostate reads Probability normalization implies that With Eqs. (9), (11) and (13) the microscopic master equation (3) can be exactly coarse-grained as follows, with the mesoscopic transition rates W N N (λ t ) = Ω N ,N w N N (λ t ). The quantity Ω N ,N takes into account that only those microstates α N and α N contribute to the sum in Eq. (14) which are connected to each other. This amounts to determine how many microstates α belong to the mesostate N under the constraint that they are connected to microstates α belonging to the mesostate N . The combinatorial problem is readily solved by noting that the occupation number that is decremented during the transition corresponds to the wanted quantity, i.e.
where N i + 1 is understood as (N i + 1) mod q. It is easy to verify that the stochastic property of the transition rate matrix is preserved by the coarse-graining, N W N N (λ t ) = 0. The mesoscopic transition rates are still consisting of multiple contributions due to the different heat reservoirs, that separately preserve the microscopic local detailed balance relation (5) at the mesoscopic level, with the notation f α,α in Eq. (5). Here, we introduced the free energy of a mesostate and used the Boltzmann entropy along with the relation which can be seen by using Eqs. (10) and (15). If the transition rates are kept constant, λ t = λ, the dynamics will reach a unique stationary state, ∂ t P s N (λ) = 0. If furthermore all heat reservoirs have the same inverse temperature, β (ν) = β ∀ν, and the nonconservative forces vanish, f (ν) = 0 ∀ν, the stationary distribution coincides with the equilibrium one which satisfies the mesoscopic detailed balance condition, and because of Eq. (17) assumes the canonical form, with the mesoscopic equilibrium free energy The marginalization of the equienergetic microstates significantly reduces the complexity of the system since the mesoscopic state space asymptotically grows like a power law, as opposed to the exponential growth of the microscopic state space. Since it will be useful further below, we remark that for a stationary mesoscopic distribution, all microstates that belong to the respective mesostates are equiprobable. This can be seen by first noting that in the stationary state, the microscopic master equation (3) reduces to 0 = j w ij p j . Since the microscopic transition rates (4) do not depend on the individual microstate α N belonging to a given mesostate N , it follows that the microscopic probability does not either in the stationary state so that A more formal proof is deferred to appendix A. We demonstrated that for thermodynamically consistent and discrete identical systems with all-to-all interactions there is an exact coarse-graining of the microscopic stochastic dynamics characterized by many-body states towards a mesoscopic stochastic dynamics that is fully characterized by the global occupation of the different unit states. It is however a priori not obvious that the thermodynamic structures built on top of these Markov process using stochastic thermodynamics are equivalent. This issue is investigated in the following section.

A. Trajectory Definitions
After having established the stochastic dynamics at microscopic and mesoscopic scales, the following is devoted to formulating the stochastic thermodynamic quantities across these scales. To this end, we first introduce the fluctuating quantities at the level of a single trajectory. Generically, a trajectory is denoted by m (τ ) (t). This notation corresponds to the specification of the actual state in the time interval under consideration, m (τ ) (t), t ∈ [t 0 , t f ]. Here, τ is a parametrization of the trajectory specifying the initial state m (τ ) (t 0 ) = α 0 , the subsequent jumps from α j−1 to α j as well as the heat reservoir ν j involved at the instances of time, t = τ j , j = 1, . . . M , and the final state, m (τ ) (t f ) = α M , where M is the total number of jumps. More explicitly, we write and refer to Fig. 1 a) for an illustrative example of such a stochastic trajectory. In the following, we will use lower scripts to label trajectory-dependent quantities in microscopic representation and write o[m (τ ) , t] for the value the observable o takes at time t for the trajectory m (τ ) . We define the energy associated with the trajectory at time t to be given by the energy of the particular microstate α the system is in for the trajectory under consideration, i.e. where the Kronecker delta δ α,m (τ ) (t) selects the state α in which the trajectory is at the time under consideration. The stochastic energy is a state function, as indicated by the notation ∆e, and its time-derivative [59] can be decomposed as follows, with the stochastic heat and work currentṡ where we introduced the notation ∇ λt = ∂ λt , ∂ λ t andẋ| m (τ ) (t) which corresponds to the instantaneous and smooth changes of x along the horizontal segments of the trajectory m (τ )(t) in Fig. 1a). It will be proven instrumental to split the fluctuating work current into the contributionẇ λ [m (τ ) , t] from the nonautonomous driving and the dissipative contribution L ν=1ẇ (ν) , t] due to the nonconservative forces. It is noteworthy that Eq. (29) is the stochastic first law and ensures energy conservation at the trajectory level [60]. As an illustrative example, Fig. 1 b) shows the time-integrated stochastic first law for the corresponding trajectory in a).
Next, the stochastic system entropy is defined as follows [61] s and is therefore also a state-function, where we set k B ≡ 1. Its time-derivative can be split into the stochastic entropy floẇ and the stochastic entropy production ratė We note that Eq. (34) corresponds to the entropy balance at the trajectory level. It will prove useful to also consider the time-integrated stochastic first law with the time-integrated fluctuating energy current and the fluctuating heat and work Using Eqs. (36) and (39), the entropy production can be written as follows B. Generating Function Techniques

Microscopic Description
In the preceding section we introduced in detail all the relevant fluctuating thermodynamic quantities. We now present techniques in order to compute the statistics and features of these quantities as they will also prove useful to determine if the thermodynamics is invariant under the dynamically exact coarse-graining in Eq. (14). To this end, we consider the microscopic generating function related to the change δo[m (τ ) , t] of the fluctuating microscopic observable o along a trajectory m (τ ) conditioned to be in a microstate α at time t which is defined as where · α denotes an ensemble average over all trajectories that are in the microstate α at time t and γ o is the counting field (also bias). It thus holds that g(γ o , t) = α g α (γ o , t). The microscopic generating function can also be expressed as follows where p(δo, t) is the probability to observe a change δo in the microscopic observable o until time t. The different moments of the microscopic observable δo are obtained via the associated microscopic generating function as follows, The equation of motion for the microscopic generating function has the form of a biased microscopic master equation [62], where w αα (γ o , λ t ) is the microscopic biased generator. The notation o α and o αα refers to the value of the stochastic microscopic observable in microstate α and its change during a transition from state α to α while the system exchanges energy with the reservoir ν, respectively.
, 0], the ensemble average over all trajectories in Eq. (42) reduces to an ensemble average with respect to the initial microstates of the trajectories only. Consequently, the microscopic generating function associated with any state function has the simple closed form Using Eqs. (27) and (32), we have for the microscopic generating functions associated with the stochastic state-like observables energy and entropy, Moreover, substituting Eqs. (30), (31), (35) and (36) into Eq. (45), we obtain for the microscopic generating functions associated with the currents

Mesoscopic Description
We rewrite the microscopic generating function (42) as follows and define the mesososcopic generating function where · α N and · N denote ensemble averages over all trajectories that are in the microstate α belonging to a given mesostate N and over all those that are in mesostate N at time t, respectively.
Thus, the microscopic generating function for the energy (47) in mesoscopic representation reads and from the microscopic equation of motion for the generating function (45) we get with the mesoscopic biased generator for O = E, Q, W, S e . More explicitly, Eqs. (49), (50) and (51) can be rewritten in mesoscopic representation as follows It is easy to verify that for O = E, Q, W, S e and o = e, q, w, s e . Thus, we find that the statistics of the stochastic first law in microscopic representation (29) is invariant under coarse-graining. Conversely, the stochastic system entropy (32) and stochastic entropy production rate (36) are functions of the microscopic ensemble probability. The corresponding equation for the mesoscopic generating function (54) would, in general, not be closed and the stochastic entropy balance in microscopic representation (34) is, in general, not invariant under the coarse-graining. Though, there are two generic cases for which an exact coarse-graining is possible.
First, for the choice of a microscopic initial condition, p sp α (0) = P sp N (0)/Ω N , where all microstates are uniformly distributed inside the respective mesostates according to Eq. (25). The local equilibrium is preserved at all times since the Hamiltonian (2) and thus the microscopic transition rates (4) do not discriminate between the equienergetic microstates inside the mesostate. For such an initial condition, the mesoscopic generating functions associated with the system entropy and entropy production rate read, respectively Secondly, according to Eq. (25), if the mesoscopic system is in a stationary state the microscopic probabilities inside a mesostate are also stationary and thus uniformly distributed regardless of any possible nonuniform distribution at intial times. Consequently, the stationary mesoscopic generating functions associated with the system entropy and entropy production rate become, respectively . Hence we conclude that the statistics of the stochastic entropy balance (34) is invariant under the coarse-graining, if one considers initial conditions which are uniform within each mesostate, or for systems in stationary states. In fact, Eqs. (25), (62) and (63) represent a potential strategy to infer the entropy fluctuations in the mesoscopic state space at finite-time: Before starting the actual measurement, the non-autonomous driving is switched off and the system is reaching a unique stationary state. The system can now be non-autonomously driven out of its steady state during the measurement and the stochastic entropies can be calculated at finite time in the mesoscopic representation via Eqs. (62) and (63).
Comparing Eqs. (47), (49) and (50) with Eqs. (56), (59) and (60), we note that the evolution of the generating functions associated with the first-law observables, that is energy, heat and work, have the same form in microscopic and mesoscopic representation. In contrast, the mesoscopic generating functions associated with the entropies do not have the same form as the microscopic ones but also contain the internal entropy S int . This is due to the coarse-grained degrees of freedom that give rise to Boltzmann entropies (19) assigned to the mesostates. Physically, the conditions for the invariance of the stochastic entropy balance [Eqs. (62) and (63) or (64) and (65)] can be understood as follows. If the microscopic degrees of freedom inside the mesostates are not equiprobable, there are transient microscopic currents that can not be grasped at the mesoscopic level and which only vanish identically once the uniform probability distributions inside the mesostates are achieved.
So far, we have established two descriptions of the stochastic thermodynamics at the microscopic and mesoscopic level. These two formulations are equivalent for the stochastic first law. In case of the stochastic entropy balance, the microscopic and mesoscopic thermodynamics coincide under the condition that the microstates inside each mesostate are equiprobable. The thermodynamics consistency at each level is ensured by the respective local detailed balance conditions in Eqs. (5) and (17). Alternatively, the thermodynamic consistency is also encoded by the so-called detailed fluctuation theorem for the stochastic entropy production. In the following, we will discuss this symmetry of the fluctuations of the entropy production as it will be of importance further below.

C. Detailed Fluctuation Theorems Across Scales
Let us consider a forward process that starts from a state that is at equilibrium with respect to the reference reservoir ν = 1, The system then evolves under the driven microscopic Markov process according to the forward protocol λ t , t ∈ [0, t].
For the backward process, indicated by the notation "˜", the system is initially prepared in the final equilibrium state of the forward process and subsequently evolves under the time-reversed driven microscopic Markov process according to the backward protocolλ t = λ t−t , t ∈ [0, t]. Then, the following microscopic finite-time detailed fluctuation theorem ensues [57,63] ln p β (1) where ∆a eq 1 = a eq 1 (λ t ) − a eq 1 (λ 0 ) denotes the change in global microscopic equilibrium free energy with respect to the reservoir ν = 1 along the forward process that only depends on the initial and final value of the driving protocol and thus does not fluctuate.
In fact, this microscopic finite-time detailed fluctuation theorem also holds for the joint probability distribution, where we write the time-integrated microscopic autonomous work currents as {δj f ) and the time-integrated microscopic energy currents as {δj is the probability to observe a microscopic nonautonomous work β (1) δw λ , the microscopic time-integrated autonomous work currents {δj  (69) can also be derived via the following symmetry of the associated microscopic generating function as demonstrated in appendix B. Analogously, we can define the forward and backward process as above also in the mesoscopic state space. In this case, the equilibrum distributions for the forward and, in reversed order for the backward trajectory, read, respectively Crucially, all fluctuating quantities appearing in the microscopic detailed fluctuation theorem (69) are invariant under the dynamically exact coarse-graining (14). Consequently, the symmetry for the microscopic generating function (70) is also exhibited at the mesoscopic level, where ∆A eq 1 = A eq 1 (λ t ) − A eq 1 (λ 0 ). Moreover, for brevity we introduced the notation and the mesoscopic time-integrated autonomous work currents {δJ f ) as well as the mesoscopic time-integrated energy currents {δJ (L) . Thus, the detailed fluctuation theorem (75) also holds at the mesoscopic level, where P β (1) E } is the probability to observe a mesoscopic nonautonomous work β (1) δW λ , the mesoscopic time-integrated autonomous work currents {δJ Having stated the various detailed fluctuation theorems across scales, we now proceed to show that the latter are relations for the entropy production of the forward processes including the relaxation from the nonequilibrium state at time t towards the final equilibrium state of the forward process (67) which coincides with the initial equilibrium state of the backward process. First, we note that initial state (66) can be prepared by disconnecting all other heat reservoirs, fixing the protocol at value λ 0 and letting the system relax. At time t = 0, all other heat reservoirs are simultaneously connected to the system and both the nonconservative f (ν) and the nonautonomous driving is switched on. As a result, the system evolves under the driven microscopic Markov process according to the forward protocol λ t , t ∈ [0, t] towards a nonequilibrium state p neq α (t). During that evolution heat δq , t] is exchanged between the system and the reservoirs ν. There is furthermore autonomous δw , t] work done on or by the system as well as nonautonomous work δw λ [m (τ ) , t] performed on the system by the external driving to change its energy landscape e α (λ t ). At time t, all heat reservoirs but the reference one ν = 1 are disconnected, the driving parameter is kept constant at its final value λ t and the nonconservative force f is switched off such that the system relaxes into the equilibrium state (67). The preparation of the starting and ending distribution of the backward process is analogous. The forward and backward process are illustrated in Fig. 2.  (67), the fluctuating entropy production (41) along the forward process can be rewritten as follows which is exactly the r.h.s of the detailed fluctuation theorem (68). Thus, the detailed fluctuation theorem (68) is a symmetry relation for the entropy production not only from initial time 0 to the time of the final protocol value t but including also the relaxation contribution from time t until the final equilibrium distribution is attained. Though, it is a finite-time relation since all fluctuating quantities in the entropy production of the forward process (76) stop evolving at time t and do thus not contribute of the statistics of the following relaxation process. We want to stress that the existence of the detailed fluctuation theorem for the entropy production across scales (69), (75) ensures that the thermodynamics formulated at each of these levels is consistent. We will make use of this result further below when we formulate the fluctuations at the macroscopic level, that is fluctuations that scale exponentially with the system size N .

D. Microscopic And Mesoscopic First And Second Law
Before turning to the macroscopic limit, for completeness, we want to formulate the thermodynamics at the ensemble level on microscopic and mesoscopic scales and hereby, because of their importance, focus on the laws of thermodynamics. Using Eq. (44) and Eqs. (47) - (50) or Eqs. (56) -(60), we arrive at the microscopic or mesoscopic first law of thermodynamics, respectively, with the average internal energy that is equivalent at microscopic and mesoscopic scale, and with the equivalent microscopic and mesoscopic heat currentṡ as well as the equivalent microscopic and mesoscopic work currentṡ Next, with Eq. (62) or (64), respectively, we find the equivalence of the average system entropy at microscopic and mesoscopic scale, Using furthermore Eq. (63) or (65), respectively, the microscopic and mesoscopic second law of thermodynamics readṡ

A. Macroscopic Fluctuations
Thus far, we have established two equivalent representations of the stochastic dynamics above, the microscopic and mesoscopic representation. We furthermore identified the conditions under which the thermodynamics at these levels coincide. In this section, the question of how to infer the fluctuations in the macroscopic limit, N → ∞, will be addressed. To shed light on this question, we will employ the Martin-Siggia-Rose formalism [65,66] which equivalently represents the Markovian jump process via a path integral. As will be demonstrated in the following, this path-integral formalism allows to establish a fluctuating description valid at macroscopic scales in the large deviation sense [67], that is for fluctuations that scale exponentially with the number of units N .
For better readability, we omit a detailed presentation of the elementary concepts underlying the construction of the path integral and refer to Refs. [47,56,68,69] where this formalism has been used in a thermodynamic context. The mesoscopic generating function G(γ O , t) associated with a mesoscopic stochastic observable O[M (τ ) , t] within the path integral representation generically reads where D[X] denotes the path-integral measure for the function X. The quantity π is the conjugated field and can be physically interpreted as the instantaneous counting field for variations in the mesostates dN . Moreover, the biased action functional L γ O [N , π] consists of the kinetic term −πṄ , the biased Hamiltonian that accounts for the current-like contributions to G(γ O , t), of a contribution due to the nonautonomous driving of a state-like contribution and of the initial condition ln P N (0). The quantities O To proceed, we notice that the microscopic stochastic entropy production with bounding Gibbs states (76) is invariant under the dynamically exact coarse-graining (14). Inserting the latter into the generic path-integral representation for a mesoscopic generating function (85), we get where we used the shorthand notation from Eq. (74). We rescale the size-extensive state variables to express them in terms of the size-intensive density n ≡ N /N . In the macroscopic limit, N → ∞, there is a single trajectory that carries all the weight of all possible paths contributing to the path integral (85). This trajectory maximizes the size-intensive action functional, max L γ [n, π] = L γ [n * , π * ], and its coordinates are therefore determined as follows where n * ≡ lim N →∞ N * /N is the continuous mean-field density. We consequently obtain via Eq. (85) the size-scaled cumulant generating function where Y denotes a vector of fields that count the size-extensive observables appearing in γ. Crucially, the size-scaled cumulant generating function associated with the entropy production (76) satisfies a symmetry that is formally equivalent to the one the mesoscopic generating function (73) exhibits, as demonstrated in appendix C. Explicitly, we have The finding of the macroscopic symmetry (92) is nontrivial since it is mathematically not obvious that the symmetry at microscopic (70) and mesoscopic (73) scales is also asymptotically preserved at macroscopic scales in spite of discarding subextensive contributions to the current statistics.
The last equation immediately stipulates the existence of a finite-time detailed fluctuation theorem in the spirit of Eq. (75) that asymptotically holds in the macroscopic limit, f /N and δE (ν) = lim N →∞ δE (ν) /N are the size-intensive equilibrium free-energy with respect to the reference reservoir ν = 1, the size-intensive nonautonomous work current, the size-intensive autonomous work current and the size-intensive energy current, respectively, in the macroscopic limit.
The existence of the finite-time detailed fluctuation theorem (93) is an important result as it ensures the thermodynamic consistency of the path-integral approach at macroscopic scales, i.e. for fluctuations that are extensive in and thus scale exponentially with the system size N .

Mean-Field Dynamics
We proceed by formulating the dynamics and thermodynamics in the macroscopic mean-field limit, where the system behaves deterministically. First we note that for an unbiased dynamics, γ = 0, that the extremal values of the auxiliary field are π * = 0. Thus, the action functional (89) needs only to be maximized with respect to the density resulting into the following Hamiltonian equations of motion The Hamiltonian equations of motion correspond to the mean-field equation governing the deterministic dynamics of the most likely occupation (mean-field) density and read explicitly, with the mean-field transition rate matrix that is stochastic, i k ij (λ t ) = 0, and whose contributions corresponding to the different heat reservoirs obey the mean-field local detailed balance ln k We note that because of probability conservation the nonlinear mean-field equation is q − 1 dimensional.

Mean-Field First And Second Law
Analogously to Sec. III D, we now want to formulate the first and second law in the macroscopic mean-field limit. Following a similar procedure as for the derivation of Eq. (95), we obtain from Eqs. (85) and (86) for the mean-field energy, whose time-derivative constitutes the first law in the macroscopic limit, with the mean-field heat and work current, A closer inspection of Eqs. (81) and (82) reveals that in the deterministic macroscopic limit the stochastic (Shannon) part of the mesoscopic system entropy vanishes and only the internal entropy of the mesostates (19) remains finite. Using the Stirling approximation we find with Eqs. (10) and (19) that the total internal entropy can be rewritten as follows, where O(ln N ) gives the order of magnitude of the error made by the approximation. As a result, using Eqs. (81) or (82) we obtain for the macroscopic entropy, respectively, where we used that (ln N )/N → 0 as N → ∞. The entropy in deterministic many-body systems therefore originates from the Boltzmann entropies related to the internal structure of the mesostates. Remarkably, the deterministic macroscopic entropy takes the form of a Shannon entropy for the mean-field density. Next, Equations (83) and (104) or Eqs. (84) and (105) stipulate for the second law in the macroscopic limit, respectively,Ṡ Hence the microscopic and mesoscopic observables in Eqs. (77)-(84) converge to the corresponding macroscopic ones in Eq. (98)-(107) if the macroscopic limit is taken, where we recall that the mesoscopic representations for O = S, Σ are only valid if the microstates inside each mesostate are equiprobable. This constitutes our main result: For thermodynamically consistent and discrete many-body systems with all-to-all interactions there is an exact coarse-graining (14) of the microscopic stochastic dynamics towards a mesoscopic one that is fully characterized by the system occupation. In the macroscopic limit, N → ∞, the stochastic dynamics asymptotically converges to a deterministic and nonlinear macroscopic (mean-field) master equation (95). Hence the stochastic dynamics can be equivalently represented across microscopic and mesoscopic scales and asymptotically on macroscopic scales as N → ∞. Furthermore, the thermodynamics can be equivalently formulated at microscopic and mesoscopic scales if the microstates inside each mesostate are equiprobable (25). The thermodynamic consistency at each of the two levels is encoded in the respective detailed fluctuation theorem, see Eqs. (69) and (75). Using a path-integral representation of the stochastic (thermo)dynamicsà la Martin-Siggia Rose, the fluctuations which scale exponentially with the system size also satisfy a detailed fluctuation theorem (93) and are therefore also thermodynamically consistent.

V. EXAMPLE
To illustrate the methodology developed in the preceding Sec. IV we consider a semi-analytically solvable autonomous Ising model which exhibits a nonequilibrium phase transition, thus representing a suitable model to demonstrate the utility of the aforementioned methods. To this end, let us consider N → ∞ spins with flat energy landscapes, 1 = 2 , that globally interact via a pair potential u/N if they occupy the same spin state i = 1, 2. The system is in contact with two heat reservoirs at different inverse temperatures β h and β c with β h < β c . According to Eq. (95), the mean-field dynamics is governed by the following nonlinear rate equation with the mean-field transition rates which we assume to be of Arrhenius form with the constant kinetic prefactor Γ that sets the time-scale of the Markov jump process. We note that the mean-field dynamics (109) is effectively a one-dimensional equation since we have n 2 = 1 − n 1 because the number of spins is conserved. We can immediately read off the stationary solution n s i = 1/2, i = 1, 2 for Eq. (109). The stability of this symmetric fixed point is encoded in the spectrum of the linearized Jacobian, A ij ≡ [∂(∂ t n i )/∂n j ]| ni,j =1/2 , which can be readily determined as follows The zero eigenvalue λ 1 reflects that the rank of the Jacobian is smaller than its dimension due to the constraint i ∂ t n i = 0. More strikingly, the second eigenvalue λ 2 changes its sign for attractive interactions, u < 0, at the critical temperatures indicative of a supercritical pitchfork bifurcation that destabilizes the symmetric fixed point into two asymmetric fixed points as can be seen in Fig. 3. This density plot depicts the stationary solution n s 1 as a function of all physical initial conditions n 1 (0) and for different cold temperatures β (c) while β (h) ≡ 1 and u = −1 are kept constant. As can be observed, the symmetric fixed point is stable for β (c) < β (c) c = 3. In contrast, for lower temperatures β (c) > β (c) c = 3 the symmetric fixed point is unstable and the system dynamics goes to one of the two asymmetric stable fixed points depending on the basin of attraction in which the initial condition lies. These two stable fixed points are related to each other via permutations of their coordinates, in agreement with the invariance of the mean-field Eq. (109) under a permutation operation. The phenomenology observed in Fig. 3 can be physically seen as follows. In the high-temperature limit the system behaves entropically, thus occupying the symmetric fixed point. Conversely, in the low-temperature limit the system behaves energetically, thus exhibiting two asymmetric fixed points that converge to the two energy ground states, that is n 1 = 1, n 2 = 0 and n 1 = 0, n 2 = 1, as β → ∞. For isothermal systems, Eq. (112) implies the critical point β c = −2/u. This is in agreement with the q-dependent universal critical temperature, β c (q) = −q/u for isothermal and all-to-all interacting q-state clock models derived in Ref. [20]. We add that the isothermal system displays a first-order equilibrium phase transition.
We now return to the non-isothermal case and consider the fluctuating quantity in Eq. (76) that for the autonomous Ising model simplifies to According to Eq. (75), our model system therefore satisfies a finite-time detailed fluctuation theorem for the timeintegrated energy current. Using the path-integral formalism introduced in Sec. IV, we however observe that analytical progress is difficult at finite time as it would require to solve the full extremization problem (90) which is analytically not possible. Instead, we therefore resort to the stationary case which considerable simplifies the problem of finding the dominant trajectory among all paths contributing to the path integral. The biased Hamiltonian (86) in the pathintegral formulation of the generating function (89) associated with the stochastic observable in the last equation reads where we added the Lagrangian multipliers λ n and λ π to enforce the spin conservation, n 1 +n 2 −1 = 0 and π 1 +π 2 = 0. The extremal value for π 1 can be solved analytically, and the extremal value n * 1 is subsequently determined numerically. In the t → ∞ limit, the boundary terms in the action functional become negligible so that the time-and size-scaled cumulant generating function is asymptotically equal to the biased Hamiltonian evaluated at the extremal values n * and π * , The scaled cumulant generating function is plotted in Fig. 4a). We choose the values β (h) = 3, β (c) = 5, u = −1 corresponding to the phase where the mean-field dynamics exhibits two asymmetric stable and a symmetric unstable fixed point. Similarly, we observe two asymmetric γ-dependent fixed points n * 1 (γ) whose coordinates are related to each other via a permutation as well as a symmetric fixed point at n * 1 = 1/2. The regime around 1/2 corresponds to the symmetric fixed point and thus to a null observable (113). Next, we note that the curve is symmetric with respect to the value γ = 1/2, thus implying that the scaled cumulant generating function asymptotically satisfies the symmetry relation which in turn stipulates the existence of a macroscopic steady-state detailed fluctuation theorem for the time-integrated energy current The existence of the steady-state fluctuation theorem is by no means obvious, here. In general, the implicit assumption underlying steady-state fluctuation theorems is that the contribution of the boundary terms related to the initial and final state of each trajectory are subextensive in time and thus negligible in the infinite-time limit. There are however situations where this may not be true, e.g. in bistable systems for starting distributions of the forward and backward . The parameters are chosen as β (h) = 3, β (c) = 5, u = −1 so that for γ (c) E = 0 the stationary mean-field system is in its energetic phase which has two asymmetric stable fixed points and a symmetric unstable one.
process that are located in the different basins of attraction. Though, in this model the two γ (c) E -dependent fixed points are related to each other via permutation of their coordinates and the statistics of the corresponding stationary states are thus identical. Figure 4b) shows the rate function Φ δJ (c) E associated with the scaled cumulant generating function G s γ (c) E ) in a). The rate function is defined as [67] Φ δJ and is related to its corresponding scaled cumulant generating function via a Legendre-Fenchel transformation Φ δJ Here, P δJ (c) E is the probability to observe a change in the energy equal to δJ (c) E and sup denotes the supremum. As can be seen in Fig. 4, both the scaled cumulant generating and the rate function are convex functions and the latter has a unique minimum equal to zero.
Our thermodynamically consistent framework allows to translate the terminology of nonlinear dynamics, i.e. the supercritical pitchfork bifurcation at the critical temperature (112), into the language of nonequilibrium statistical mechanics, i.e. a nonequilibrium phase transition at the same critical temperature. For this purpose, we prepare the system in its critical state by setting β (h) = 1, β (c) = 3, u = −1. Fig. 5 depicts in a) the scaled cumulant generating function (117) with the system being in its critical state. The scaled cumulant generating function exhibits a kink at γ (c) E = 0 indicative of a nonequilibrium phase transition. Owing to the symmetry (118), the scaled cumulant generating function has another kink at γ (c) E = 1. The non-differentiability of the generating function at γ (c) E = 0 implies that the rate function in Fig. 5 b) would be nonconvex over a finite interval. The Legendre-Fenchel transformation (121) yields not the nononvex rate function but its convex envelope Φ ce δJ (c) E . Here, the part of the convex envelope that replaces the nonconvex regime of the rate function corresponds to the flat part of the curve in the vicinity of the δJ (c) E = 0. Thus, we find that the time-integrated energy current distribution in Eq. (120) is bimodal, thus also encoding the nonequilibrium phase transition.

VI. CONCLUSION
In this work we demonstrated how to consistently build a stochastic dynamics and thermodynamics description across scales for many-body systems with all-to-all interactions: For this purpose, we considered a system of N allto-all interacting identical and classical units consisting of q states. The units undergo transitions due to several heat reservoirs and because of external forces. The microscopic stochastic dynamics characterized by many-body states can be exactly coarse-grained towards a mesoscopic one that is determined by the occupation numbers of the different unit states. Here, the all-to-all interactions give rise to equienergetic many-body states which form the mesostates. Importantly, the coarse-graining significantly reduces the complexity of the many-body system as the growth of the state space changes from an exponential to a power-law one. Employing the formalism of stochastic thermodynamics, it was proven that the stochastic first law of thermodynamics is always invariant under the dynamically exact coarsegraining. Conversely, this only holds true for the stochastic entropy balance if the microstates within each mesostate are equiprobable.
We then considered the macroscopic limit, N → ∞. To consistently determine the macroscopic fluctuations we used the Martin-Siggia-Rose formalism. We showed that the fluctuations that scale exponentially with the system size N are thermodynamically consistent as they obey a detailed fluctuation theorem. Detailed fluctuation theorems of the same form were also derived at the microscopic and mesoscopic level, hence proving thermodynamic consistency across scales. Moreover we proved via the path integral representation of the stochastic dynamics that the mesoscopic master equation asymptotically converges to a nonlinear rate equation. The methodology to determine macroscopic fluctuations was demonstrated via a semi-analytically solvable Ising model in contact with two reservoirs and exhibiting a nonequilibrium phase transition. Our work provides a powerful framework to address the thermodynamics of nonequilibrium phase transitions.
An interesting outcome of this work is that the thermodynamic description of many-body all-to-all interacting systems, when going from a microscopic to an occupation level description, assigns Boltzmann entropies (logarithms of complexion numbers) to each mesostate, despite the fact that the system is driven away from-equilibriun by multiple reservoirs and external forces. Furthermore, in the deterministic macroscopic limit, N → ∞, the ensuing entropy takes the form of a Shannon entropy for the deterministic occupation which exclusively results from these internal mesostate entropies.
We end by placing our findings in the context of the recent works on thermodynamically consistent coarse-graining. Many of these are based on time-scale separation: fast degrees of freedom reach a local stationary state over timescales much shorter than the slow dynamics and can be adiabatically eliminated. The resulting transition rates of the slow dynamics then satisfy a local detailed balance condition which carries the information about the thermodynamic potentials (energetic and/or entropic) [70][71][72] or the driving forces [73][74][75] resulting from the fast dynamics. Some coarse-grainings do not require time-scale separation and the hidden degrees of freedom have been shown to behave as work sources (pure energy no entropy) on the remaining degrees of freedom, see e.g. Refs. [72,76]. In the present work, the coarse-graining does not rely on any time-scale separation but results from the all-to-all interactions which do not discriminate energetically between the different microstates leading to the same global occupation. As a result a purely entropic contributions ensues at the occupation level. Models where such coarse grainings appeared can be found in Refs. [18,20,21]. We finally present the proof of the asymptotic symmetry in Eq. (92). First, we rescale the extensive state variables n = N /N so that N O n ≡ O(N ). Also, the discrete gradient of an observable along an edge asymptotically becomes a derivative with growing size, N (∂ ni − ∂ nj ) O n = O ij (N ). The path-integral representation of the generating function (89) then reads N (t )) .

(C1)
The crucial step of the derivation is to define physically consistent transformation rules to time-reverse the biased stochastic dynamics. Time-reversal transformations of unbiased Langevin dynamics have been investigated in Ref. [77]. For the generating function in question (C1), we define the time-reversed biased stochastic dynamics as follows t = t − t ,ñ = n,λ t = λ t−t , π = −π + N β (1) ∇ n A (1) (n) = −π + N β (1) ∇ n E(n) − N ∇ n S int (n),γ = 1 − γ, while reusing the shorthand notation from Eq. (74) and, for better readability, rewrite the state functions in a more compact way as follows O n (λ t ) → O(n). The definitions of the time-reversed physical quantities in the first line are trivial. Less obvious is the transformation rule of the auxiliary field π. This transformation rule amounts to inverting the directions of the edges corresponding to a reversion of the Markov dynamics: The change of the sign in front of π can be seen by noting that the latter is a counting field for variations in the state variables dN . Moreover, the affinity along an edge is inverted by the free energy shift. We proceed by demonstrating that the above transformation, up to a non-fluctuating quantity, indeed leaves the generating function invariant. For better readability, we will split the action functional (C1) into two parts and investigate how they transform under the time reversal in Eq. (C2). First, the invariance of the biased Hamiltonian under this time-reversal transformation can be seen as follows, exp πj(t ) − πi(t ) − N β (1) (∂n j − ∂n i )A (1) (n(t )) + 1 − γ (ν) Furthermore, we find for the sum of the kinetic and non-autonomous driving terms together with the initial condition under time-reversal, t 0 dt N β (1) ∇ n A (1) (n(t )) − π(t ) ·Ṅ (t ) + N (1 − γ Λ ) β (1)λ t · ∇ λ t E n (λ t ) + ln P eq Nt (λ t ) = = t 0 dt N β (1) d t A (1) (n(t ))−λ t · ∇ λ t A (1) (n(t )) −π(t )·Ṅ (t )+N (1−γ Λ ) β (1)λ t · ∇ λ t E n (λ t ) +ln P eq Nt (λ t ) = − t 0 dt {π(t ) ·Ṅ (t ) + N γ Λ β (1)λ t · ∇ λ t E n (λ t ) + N β (1) A (1) (n(t)) − A (1) (n(0)) + ln P eq Nt (λ t ) Collecting results, we thus find that the size-intensive action functional is invariant under the time reversal (C2) up to a non-fluctuating term corresponding to the change in the size-intensive part of the equilibrium free energy, i.e.
L γ [n, π] =Lγ[n, π] − β (1) ∆A eq 1 (λ). (C5) In the macroscopic limit, the scaled cumulant generating function is equal to the extremal action functional, cf. Eq. (91). Moreover, the action functional contains the initial condition of the trajectories so that its extremization does not give rise to additional boundary terms, Hence the invariance of the action functional is preserved in the macroscopic limit that in turn stipulates the following symmetry for the scaled-cumulant generating function, which is exactly Eq. (92).