Thermodynamics of interacting systems: the role of the topology and collective effects

We will study a class of system composed of interacting unicyclic machines placed in contact with a hot and cold thermal baths subjected to a non-conservative driving worksource. Despite their simplicity, these models showcase an intricate array of phenomena, including pump and heat engine regimes as well as a discontinuous phase transition. We will look at three distinctive topologies: a minimal and beyond minimal (homogeneous and heterogeneous interaction structures). The former case is represented by stark different networks ("all-to-all"interactions and only a central interacting to its neighbors) and present exact solutions, whereas homogeneous and heterogeneous structures have been analyzed by numerical simulations. We find that the topology plays a major role on the thermodynamic performance for smaller values of individual energies, in part due to the presence of first-order phase-transitions.Contrariwise, the topology becomes less important as individual energies increases and results are well-described by a system with all-to-all interactions.

This paper aims to partially fill this gap, by studying the influence of topology of a simple class of system, composed of interacting units, referred to here as a collection of nanomachines, placed in contact with a hot and cold ther-mal baths subjected to a non-conservative driving worksource.The approach to be considered here is akin to the commonly referred to as "lattice-gas" models in the realm of equilibrium statistical mechanics.They have a longstanding importance in the context of collective effects and serving as the cornerstone for numerous theoretical, experimental, and technological breakthroughs, encompassing the ferromagnetism, liquid phases, the topology effect, the fluctuation-driven generation of new phases and others, highlighting that distinct systems have been described/characterized via Hamiltonian of the fundamental models (e.g., the Ising, Potts, XY and Heisenberg models).The all-to-all version for our model has been investigated previously [53,54,58] for finite and infinite number of interacting units, in which the cooperative effect gives rise to a rich behavior, including the enhancement of the power and efficiency at optimal interactions, the existence of distinct operation models and a discontinuous phase transition.Furthermore, systems with all-to-all interaction can be solved analytically, which makes them easier to analyze.There are, however, also many systems where systems only interact locally (nearest-neighbor like).
In this paper, we present a detailed study on the influence of the topology in aforementioned class of interacting units.We will focus on two distinctive approaches: minimal models, which can be treated analytically, and more complicated systems, where we will focus on numerical analysis.In the former class, we will focus on systems with all-to-all interaction and a one-to-all interaction (also known as stargraph), in which a single central spin is interacting with all other units.After that, we go beyond the minimal models by considering the influence of homogeneous and heterogeneous interaction topologies.We will show that, for small values of individual energies β ν ϵ of each occupied unit, the lattice topology can have a significant influence on the system performance, in which the increase of interaction V among units can give rise to a discontinuous phase transition.Conversely, as β ν ϵ increases, the phase transition is absent and the topology plays no crucial role and the models seem mutually similar.In this case one can get approximate expressions for thermodynamic quantities through a phenomenological two-state model.
The paper is structured as follows: in Section II, we introduce the model and its thermodynamics.Section III describes the lattice topologies which will be analysed.In Section IV, the aforementioned minimal models are studied.In Section V, the more complicated topologies are studied.Conclusions are drawn in Section VI.

II. MODEL AND THERMODYNAMICS
We are assuming a system composed of N interacting twostate nanomachines.The two states of each individual machine are denoted as σ i = 0(1) according to whether it occupies the lower(upper) state with energy 0(ϵ).We will consider arXiv:2308.02255v2[cond-mat.stat-mech]17 Oct 2023 that the system is in contact with two thermal baths at different temperatures.Furthermore, there will be a non-conservative force (described below) that extracts work from the system, in this way creating a thermal engine.The state of the full system is then described by σ ≡ (σ 1 , σ 2 , ..., σ i , ..., σ N ), where σ i describes the state of the i'th machine.Throughout this paper, we shall restrict our analysis on transitions between configurations σ and σ i differ by the state of one machine, namely that of unit i.In this case, the time evolution of probability p(σ, t) satisfies a master equation, where σ i ≡ (σ 1 , ..., 1 − σ i , ..., σ N ) and index ν = 1(2) accounts for transitions induced by the cold (hot) thermal bath.The transition rate due to the contact with the ν-th thermal bath are assumed to be of Arrhenius form where ∆E i (σ) is the difference of energy between states σ i and σ and Γe −β ν E a /2 accounts to the coupling between the QD and thermal bath, expressed in terms of the activation energy E a and β ν = 1/T ν , [hereafter we shall adopt k B = 1].As stated previously, the interaction among units will depend on the lattice topology, whose energy of system is given by the generic expression: where n = N i=1 σ i denotes the total number of units in the state of energy ϵ, V is the interaction strength (δ being the Kronecker delta) and ⟨k⟩ is the average number of neighbours to which each unit is connected.Eq. (3) has been inspired by earlier studies about interacting system, in which a similar type of interaction is consider to describe the interaction between nanomachines in distinct states [53,54].Also, this interaction shares some similarities with recent papers [59,60] in which the tunneling between two quantum-dots is investigated via the inclusion of a similar term.We will both look at topologies where k i , the number of nearest neighbors of the unit i is independent of i, ⟨k⟩ = k i (all-to-all interactions and homogeneous systems) and cases where k i depends on i (stargraph and heterogeneous systems).One of these earlier studies also used similar types of work sources: we consider the worksource given by F ν with F ν = (−1) ν (1 − 2σ i )F, in such a way that transitions 0 → 1 (1 → 0) are favored according to whether the system is placed in contact with the cold (hot) thermal baths.This type of interaction can also be mapped on other types of systems such as kinesin [61,62], photo-acids [63] and ATP-driven chaperones [64].
From Eq. (1) together transition rates given by Eq. ( 2), the time evolution of mean density ρ = ⟨σ i ⟩ and mean energy ⟨E(σ)⟩ = σ E(σ)p(σ, t) obey the following expressions: ρ = ⟨(1 − 2σ i )(ω (1)  i (σ and respectively, where P and ⟨ Qν ⟩ denote the mean power and the heat exchanged with the ν-th thermal bath and are given by [53,57]: and the standard stochastic thermodynamics expressions for power and heat respectively [3].Throughout this paper, we will assume that the system has relaxed to a steady state, p(σ, t) → p st (σ), i.e., P + ⟨ Q1 ⟩ + ⟨ Q2 ⟩ = 0.In this case, one can also write the entropy production as with One can verify from Eq. ( 2) that the entropy production, Eq. ( 8) assumes the classical form Under the correct choice of parameters, an amount of heat extracted from the hot bath ⟨ Q2 ⟩ > 0 can be partially converted into power output P < 0 and the system can be used as a heat engine.The efficiency is then defined as η = −P/⟨ Q2 ⟩, which satisfy the classical relation η ≤ η c = 1 − β 2 /β 1 .

III. LATTICE TOPOLOGIES
As stated in the previous section, we intend to study the differences in thermodynamic quantities between topologies.We will focus on two classes of systems: Minimal structures, namely stargraph and all-to-all interacting systems, and beyond minimal structures, comprising homogeneous and heterogeneous systems.In stargraph systems the interactions are restricted to a central unit(hub) which interacts with its all nearest neighbor sites (leaves).
Homogeneous and heterogeneous structures present remarkably different properties and has been subject of extensive investigation.While the former case have been largely studied for addressing the main properties of graphs, the latter describes a broad class of systems, such as ecosystems, the Internet, the spreading of rumors and news, citations and others, in which the agents form heterogeneous networks and are approximately scale-free, containing few nodes (called hubs) with unusually high degree as compared to the other nodes of the network.For the homogeneous case, we shall consider those characterized by a fixed neighborhood per unit, being grouped out in two categories, including a regular arrangement (interaction between nearest neighbors) or a randomregular structure, in which all units have the same number of nearest neighbors, but they are randomly distributed.Such latter case is commonly generated through a configurational by Bollobás [65].Finally, among the distinct heterogeneous structures, we will consider the Barabasi-Albert scale-free network, being probably the most well-known example of heterogeneous networks [66].The Barabási-Albert (BA) model is based on a preferential attachment mechanism, in which the degree distribution follows a power-law with scaling exponent γ = 3 [66].

IV. MINIMAL MODELS FOR COLLECTIVE EFFECTS: ALL-TO-ALL INTERACTIONS VERSUS STARGRAPH
We will first look at the thermodynamic properties of "allto-all" and stargraph minimal topologies.There are several reasons for this.First, both of these models can in principle be solved exactly.Second, these structures can be seen as each others opposite.Third, the thermodynamic properties stargraph topologies can give some insights about heterogeneous networks (e.g.Barabasi-Albert), in a which some nodes are highly connected and most the remaining ones have few connections [31,67,68].

A. Steady-state distribution
For an all-to-all topology, the state of the system is fully characterized by the number of units in the upper state, n = N i=1 σ i .In terms of total occupation, the master equation for the all-to-all system simplifies to The steady-state distribution for p st (n) then satisfies [53] where Z is the normalization factor and ω i j = 2 ν=1 ω (ν) i j , with transition rates solely expressed in terms of n by The thermodynamic quantities can be evaluated from the probability distribution such that, expressed in terms of the probability current ).An overview of the model features in all-to-all topologies will be depicted next (see e.g.Figs. 2 and Refs.[53,54]), being strongly dependent on the interplay between individual ϵ and interaction V parameters.For β ν ϵ << 1, the increase of interaction strength V favors a full occupation of units in the upper state ρ → 1, whereas ρ exhibits a monotonous decreasing behavior upon V is raised for β ν ϵ >> 1.The crossover between above regimes yields for finite ϵ and depends on E a , β 1 /β 2 and F. Another important point to be highlighted concerns that as E a is increased and β ν ϵ is small, the interaction marks two distinct trends of the density: its decreasing behavior of prior the threshold interaction followed by sharp increase towards ρ → 1 at V = V 0 [see also e.g.Fig. 2(a)].Such behavior corresponds to a discontinuous phase transition (see e.g.Fig. 6)(a).In Sec.IV B, we shall investigate these consequences over the system performance.
It is in principle possible to calculate the exact steady-state distribution for a finite-size star-graph by diagonalising the evolution matrix.However, some insights into the dynamics can be obtained by doing appropriate approximations, as we will show now.First, we note that the state of the system can be written in terms of n and c, denoting the number of leaves and the hub in the upper state, respectively.The associated probability distribution, p(n, c, t), satisfies where and with n = 0, 1, ..., N − 1 and transition rates rewritten in the following way We assume that the hub dynamics evolves into a faster timescales than the relaxation of the surrounding leaves, in such a way that it can be assumed/treated as thermalized at a local leaf transition n ± 1.In other words transitions are such that (κ (n,1)  0,1 + κ (n,2) 0,1 )p(1|n) = (κ (n,1) 1,0 + κ (n,2) 1,0 )p(0|n) and hence, the joint probability p st (1|n) is given by: where p st (0|n) = 1 − p st (1|n).By summing Eq. ( 13) over c, together the property p(n, c) = p(c|n)p(n), one arrives at the following equation for the time evolution of probability p(n, t) where Since Eq. ( 20) is analogous to Eq. ( 9) for the all-to-all case, the probability distribution of leaves p st (n) is given by in which π i, j ≡ 2 ν=1 π (ν) i, j and Z is again the corresponding normalization factor.Thermodynamic properties can be directly evaluated from Eqs. ( 19) and ( 22), such as the system density given by ρ = N−1 n=0 [n+p st (1|n)]p st (n)/N, where the probability p h of hub to bein the upper state with individual energy ϵ reads p h = N−1 n=0 p st (1|n)p st (n).As previously, from the probability distribution, thermodynamic quantities are directly evaluated and given by where n,n+1 p st (n + 1) and p st (c, n) = p st (c|n)p st (n) was considered.Likewise, each heat component Qν is given by where J (c,ν) n+1,−1 and K (n,ν) 1 are evaluated from Eqs. ( 14) and ( 15) in the NESS.
Fig. 1 compares the evaluation of system density ρ from the exact (continuous lines) method with the approximate (symbols) method.Both curves agree remarkably well.

B. General features and heat maps for finite N
To reduce the number of parameters, we will assume that β 2 = 1 unless specified otherwise.Furthermore, we will look at ϵ = 0.1 and 1 and vary β 1 .Figs. 2 and 3 summarize the main findings about minimal models for interacting for a small system of size N = 20.
As the all-to-all, the increase of interaction strength V also changes ρ significantly for the stargraph and, consequently, affects the engine performance.While intermediate densities favor the system operation as a pump, their emptying when V FIG. 1.For the stargraph, the comparison between the exact system density ρ, obtained by diagonalising the evolution matrix (continuous lines), and the approximation Eqs. ( 19) and ( 22) (symbols).Circles and stars correspond to β 1 = 5 and β 1 = 10, respectively.Parameters: is increased changes the operation regime, from a pump to a heat engine and also increases the engine performance, whose performances are meaningfully different for ϵ = 0.1 (smaller) and 1 (larger) individual β ν ϵ's.The maximal reachable efficiency η ME is always lower than η c for finite N, as expected.
Another common behavior in both cases is the fact that large V favors a full occupation of the upper state when β ν ϵ is small (see e.g.panels (a) in Figs. 2 and 3), implying that the system operates dudly when most of units are in the upper state, whose crossover from heat to dud regime is marked by a discontinuous phase transition.Conversely, the increase of ϵ marks the absence of phase transition for a broader range of V and consequently not only extends the engine regime but also improves system performance.Despite closely dependent on parameters, both η and P exhibit similar trends as β 1 is raised for the all to all case.Although having inferior performance than the all-to-all (at least for the chosen set of parameters), the stargraph yield some striking features for smaller N (see e.g. in Figs. 3, 5 and 6), including an intermediate sets of V in which both η and P do not behave monotonously, characterized by a local and global maximum (η MP ) and minimum (P mP ), as can be seen in Figs. 3 (e)−( f ).In all cases, exact results (continuous) agree very well with numerical simulations (symbols).
A global phase-portrait is depicted in Figs. 4 and 5 for N = 20.These results are in agreement with the aforementioned and reinforce previous findings, including larger maximum efficiencies and power for all-to-all interactions than stargraph ones for small N's, but such later one presents two  distinct regions (for lower and larger V's) which the heat engine operates more efficiently.Similar findings are shown in Sec.A for P's.
In Secs.IV C and IV D, remarkable aspects about both minimal structures, including the existence of a discontinuous transition for smaller individual energies as well as its suppression as ϵ increases, shall be described.

C. Effect of system sizes and phase transitions
The first common aspect regarding the behavior of stargraph and all-to-all interaction structures is that the increase of interaction V (for smaller values of β ν ϵ) not only influences the system properties and the engine's performance but also gives rise to a phase transition characterized by a full occupancy of units in the upper state as N → ∞.However, the  behavior of finite systems provides some clues about the classification of phase transition, as described by the finite size scaling theory [69][70][71][72][73][74][75].In the present case, the existence of a crossing among curves for distinct (finite) system sizes N reveals a discontinuous phase transition [74,75], as depicted in Fig. 6.More specifically, the density curves ρ strongly depend on the system size near phase-transition V 0 (V 0 = 3.712(3) and 1.942 (2) for the all-to-all and stargraph, respectively), whose intersection among curves is consistent to a density jump for N → ∞.Such features are also manifested in the behavior of both η and P (see arrows in Figs.6), marking the coexistence between heat engine and dud regimes.The opposite scenario is verified by raising ϵ, as depicted in Fig. 7 for ϵ = 1 for both all-to-all and stargraph topologies.Unlike the behavior of ϵ = 0.5, the phase transition is absent for both structures and as a consequence, the heat engine regime is broader.
We close this section by stressing that, although discontinuous phase transition have already been reported for similar FIG. 4. For the all-to-all topology, the efficiency heat maps for various choices of ϵ' as in Fig. 3 systems [53], the existence of a phase transition in the stargraph structure is revealing and suggests that a minimal interaction structure is sufficient for introducing collective effects that are responsible for the phase transition.

D. The N → ∞ limit and phenomenological descriptions
A question which naturally arises concerns the system behavior in the thermodynamic limit N → ∞ for both all-toall and stargraph structures.The former case is rather simple and can be derived directly from transition rates, in which system behavior is described by a master equation with non linear transition rates.Since the all-to-all dynamics is fully characterized by the quantity n, the macroscopic dynamics is given by the probability of occupation p 1 , corresponding to ρ in the thermodynamic limit p 1 = lim N→∞ N i=1 ip st (i)/N (p 0 = 1 − p 1 ) [53,57] and described by the master equation that has the form ṗ1 = 2 ν=1 J (ν) 10 , such that: where transition rates ω (ν) 10 and ω (ν) 01 denote the transition the lower to the higher stateand vice versa, respectively, and are listed below: where ∆E 10 = V(1 − 2p 1 ) + ϵ.For N → ∞, expressions for the power ⟨P⟩ ≡ P/N and heat per unit ⟨ Qν ⟩ ≡ ⟨ Qν ⟩/N from Eqs. ( 11) and ( 12) read ⟨P⟩ = F(J (1) 10 − J (2) 10 ) and ( 27) respectively.p 1 is obtained by solving Eq. ( 25).As shown in Sec.IV C, small and large values of individual energy β ν ϵ mark different behaviors a N increases, the former yielding a discontinuous phase transition.Unlike the behavior of finite N, the discontinuous phase transition is featured by the existence of a hysteretic branch in which the system has a bistable behavior [53].We shall focus on ϵ = 1 which describes the behavior of large β ν ϵ's, as depicted in Fig. 7, together with a comparison with different N's.As can be seen in this figure, in both cases, maximum efficiencies η ME 's (for coupling strength V=V ME 's) and (absolute) minimum powers P mP 's (for coupling strength V = V mP 's) increase as N is raised and approaching to the N → ∞, consistent with enhancing collective effects.However, contrasting with the power, the efficiency for smaller system sizes is larger for V > V ME .This can be understood from the interplay between power and ⟨ Q2 ⟩.For V > V mP , the power mildly changes with N, whereas ⟨ Q2 ⟩/N monotonically increases with N. Likewise for V ME < V < V mP , but in this case ⟨ Q2 ⟩/N increases "faster" than ⟨P⟩.
Although Eq. ( 25) can be solved numerically for any set of parameters, its nonlinear shape makes it impossible to obtain analytical results.However, it is possible to get some insights about the system in the heat engine regime when p 1 << 1 (and p 0 is close to 1).In this case, the terms p 1 and p 0 appearing in transition rates can be neglected and treated as p 0 ≈ 1, respectively, in such a way one arrives at the followingformula: By inserting transition rates from Eq. ( 26) into Eq.( 29), we arrive at the following approximate expression for p 1 : where A ν = e − βν 2 E a .Approximate expressions for P and ⟨Q ν ⟩'s in the heat engine are promptly obtained inserting Eq. (30) into Eqs.( 27) and (28), respectively.Although they are cumbersome, they solely depend on the model parameters β 1 , β 2 , E a , ϵ, F and V.The comparison between exact and approximate results is also shown in top panels from Fig. 7 (symbols) for ϵ = 1, in which no phase transition yields (at least for limited V's).As can be seen, the agreement is very good for p 1 << 1.
We pause again to make a few comments about Eq. (36).For small values of β ν ϵ, in which a discontinuous phase transition yields at V 0 , p h jumps from 0 to 1 for V < V 0− and V > V 0+ , respectively, and p h = 1/2 precisely at V = V 0 .Second, from the hub behavior, it follows that p 1 jumps from p 1− to p 1+ , where p 1− (p 1+ ) are obtained from Eq. ( 31) evaluated at c = 0 (for V → V 0− ) and c = 1 (for V → V 0+ ), respectively.Third, the order-parameter jump is also followed by discontinuities in the behavior of thermodynamic quantities, such as ⟨P⟩ and η.They are evaluated from Eqs. ( 34) and ( 36) at c = 0 and V = V 0− to c = 1 and V = V 0+ , respectively.Fourth and last, large β ν ϵ's mark no phase transitions for limited values of V and hence the power and efficiency are evaluated at c = 0 (since p h = 0).All above findings, together the reliability of Eqs. ( 31), ( 34) and ( 36) for small and large values of β ν ϵ, are depicted in bottom panels from Figs. 6 and 7 for ϵ = 0.5 and 1, respectively.As for the all-to-all, thermodynamics quantities for finite N approach to the N → ∞ limit as N is increased.
We close this section by drawing a comparison between allto-all and stargraph performances in the heat engine regime for the same set of parameters in Figs. 6 and 7.While the former structure is more efficient, stargraph ones present larger power outputs.

V. BEYOND THE MINIMAL MODELS: HOMOGENEOUS AND HETEROGENEOUS TOPOLOGIES
In this section, we will go beyond the minimal models and look at both homogeneous and heterogeneous structures.Unlike minimal models, it is not possible to obtain analytical expressions and our analysis will focus on numerical simulations using the Gillespie method [76].Due to the existence of several parameters (β 1 , β 2 , F, ϵ, V and E a ), we shall center our analysis on F = 1, β 1 = 5 and β 2 = 1, in which results minimal models predict a marked engine regime as V is varied.Fig. 8 and 9 depict some results for homogeneous and heterogeneous structures, respectively.
Starting our analysis for ϵ = 0.1 (top panels) and homogeneous arrangements, we see [Fig.8(a) and (c)] that the system performance increases by increasing the connectivity k and there are small differences between regular and randomregular arrangements.Unlike the homogeneous case, differences between ⟨k⟩'s are particularly clear for heterogeneous structures, where the heat engine is absent for ⟨k⟩ = 4.In this case, the system only operates as a pump, similarly to the stargraph, see e.g.Fig. 3 for N = 20.On the other hand, the heat engine is present for ⟨k⟩ = 10 and 40.A possible explanation is that the former and latter cases are closer to the stargraph and the all to all structures, respectively.
The results for ϵ = 1 (bottom panels) are remarkably different.We see that the heat engine regime becomes much larger in terms of V (with ρ monotonously decreasing as V goes up).Furthermore, one can see that the influence of the lattice topology and neighborhood becomes negligible and the results become very similar to those of the all-to-all topology, revealing that the role of topology is not so important for larger β ν ϵ's.

VI. CONCLUSIONS
In this paper, we the role of topology of interactions on the performance of thermal engines.We investigated four distinct topologies for a simple setup composed of interacting unicyclic machines, each one allowed to be in two states: all-to-all, stargraph, homogeneous and heterogeneous structures.Different findings can be extracted from the present study.Interestingly, the interplay among parameters (individ- ual β ν ϵ, interaction energies V and temperatures) provides two opposite scenarios, in which the role of topology is important and less important respectively, depending on whether β ν ϵ is small or large.The former case not only shows a discontinuous phase transition as the interaction is raised, but also how the increase of neighborhood (both homogeneous and heterogeneous) increases the efficiency but in contrast its power is inferior.Since a majority fraction of are empty in the latter case, the topology of interactions plays no major role.
As a final comment, we mention some ideas for future research.It might be interesting to study the the full statistics of power and efficiency in different lattice topologies, in order to tackle the influences of fluctuations.Also, it might be interesting to compare the performance of different engine projections, such as those composed interacting units placed in contact with only one thermal bath in instead of two, in order to compare the system's performances as well a the influence of lattice topology in those cases.Finally, it shall be interest-
FIG.4.For the all-to-all topology, the efficiency heat maps for various choices of ϵ' as in Fig.3.HE, P and D denote the heat engine, pump and dud regimes, respectively.Parameters: N = 20, E a = 2, β 1 = 10, β 2 = 1.

1 FIG. 6 .
FIG. 6.The effect of system size in minimal collectively models: Left and right panels depict the behavior of density (top), η/η c (center) and ⟨P⟩ ≡ P/N (bottom) for the all-to-all and the stargraph, respectively.Arrows indicate the discontinuous phase transitions, characterized by the crossing among curves.The inset in the right bottom panel indicates the hub density for the stargraph model.Dashed lines: Results for N → ∞.Parameters: F = 1, E a = 2, ϵ = 0.5, β 1 = 10 and β 2 = 1.