Stochastic thermodynamics in many-particle systems

We study the thermodynamic properties of a microscopic model of coupled oscillators that exhibits a dynamical phase transition from a desynchronized to a synchronized phase. We consider two different configurations for the thermodynamic forces applied on the oscillators, one resembling the macroscopic power grids, and one resembling autonomous molecular motors. We characterize the input and the output power as well as the efficiency at maximum power, providing analytic expressions for such quantities near the critical coupling strength. We discuss the role of the quenched disorder in the thermodynamic force distributions and show that such a disorder may lead to an enhancement of the efficiency at maximum power.


Introduction
The stochastic thermodynamics of microscopic systems has been the subject of intense investigation in recent years [1] in an attempt to extend basic concepts of macroscopic classical thermodynamics to the microscopic realm in general, and to out-of-equilibrium microscopic systems in particular. Notably a lot of effort has been devoted to the characterization of the efficiency of microscopic devices, that can transform heat or chemical energy into mechanical work. While for the first type of devices (autonomous heat engines) the efficiency is bounded by the Carnot limit [2,3], in the case of isothermal engines the efficiency is constrained by the thermodynamic limit 1. In both cases the upper limit is reached for quasi static operation, resulting in a vanishing power output. Thus, a more relevant quantity to study is the efficiency at maximum power (EMP), that exhibits an interesting universal behaviour for different types of devices [2,3,4,5,6,7,8].
In particular the EMP in the linear regime is 1/2 of the maximal allowed value, while the behaviour beyond the linear regime depends on the details of the coupling between the energy producing and the energy consuming cycles.
Many of the theoretical studies have been directed toward the characterization of the EMP in single devices such as soft nanomachines [5], single molecular motors [6,7], devices involving single electron transport [4,9,10], or single entropy-driven motors [11]. However, an important class of microscopic devices is represented by cell molecular motors, which operate in crowded environments where their mutual interaction can become significant. For example, many kinesin motors walk on the same microtubule leading to traffic jam formation in some case [12], while several motors can pull the same cargo resulting in a strong cooperative effect [13,14,15,16,17]. Furthermore, recent studies on synthetic nanomotors have been conducted with the aim of reproducing the performance of their biological counterparts [18,19,20]. In this regard, it is important to note that it is now possible to engineer molecular spiders that exhibit directional movement and behave like robots by carrying out a sequence of predetermined actions [20]. These artificial motors might in the future be organized in teams, to optimize, for example, their transport properties and efficiency [21]. Thus it is clear that future research on thermodynamic property optimization will deal with teams of interacting motors, where the dynamical phase these motors operate in becomes relevant.
One of the most studied models of interacting particles in out-of-equilibrium physics is the exclusion process, which exhibits three distinct dynamical phases with different densities and particle currents [24]. Furthermore, the exclusion process is often used to model molecular motors moving on a lattice, see, e.g. [25,26]. We have previously investigated the issue of EMP in isothermal interacting motors, modelled as an exclusion process on a single lattice [8,22], or on a network [23]. In these studies, we found an increase of the EMP in a many-motor system with respect to the single motor case, for a suitable choice of the model parameters. Remarkably, in [8,22] we found that the enhancement of the EMP occurs in a range of parameter values compatible with the biological estimates for the molecular motor Kinesin. From those studies we concluded that after a dynamical phase transition the dynamical response of the system to an external drive can change, leading in turn to a change in the thermodynamic properties. Specifically the dependence of the delivered power on the driving thermodynamic forces may vary.
One of the main limitation that one faces when studying the thermodynamic properties of exclusion processes on a lattice, is that the intensity of interaction between the motors can only be indirectly tuned by changing the kinetic parameters and thus the density of motors on the lattice [8,22,23]. Furthermore, if one wants to study the effect of force disorder on the motor particles, one has to resort to numerical simulations, as no exact result exists for the exclusion process with heterogeneous particles.
Instead, here we consider a model of N interacting microscopic particles, where the particle-particle interaction is an explicit parameter, that can be tuned in order to drive a dynamical phase transition, from a weakly interacting -incoherent system to a strongly interacting -coherent system. This model was originally introduced by Sakaguchi in [27] as an extension of the Kuramoto model (KM) [28], to study the synchronization of a group of interacting oscillators in contact with a reservoir at constant temperature. Furthermore, the effect of quenched disorder in the thermodynamic force distribution can be taken into account within the present model. The model is introduced and discussed in section 2. The N interacting particles can be viewed as a network of energy producers and users or as a system of interacting autonomous motors under the effect of thermodynamic forces. Since the dynamical phase diagram can be obtained in terms of the particle interaction strength, temperature and force distribution, in section 4 we will discuss how to calculate the relevant thermodynamic quantities, namely the deliverer and input power, and the efficiency. We will consider two possible scenarios as far as the force distribution is concerned. In the first one, sec. 4.1, either a positive or a negative force is applied on each particle. In the second scenario, sec. 4.2, both a positive and a negative force is applied on the same particle. We will thus discuss how to optimize the delivered power for the different types of network models, hence obtaining the EMP in terms of the interaction intensity, and thus of the coherence between the particles' motion. We will finally discuss the effect of the quenched disorder in the force distribution on the thermodynamic quantities.
Interestingly, the model that we use here is a microscopic version of a model used to mimic macroscopic power grids. Indeed the dynamics of interconnected power grids can be mathematically represented by a complex network of coupled oscillators [29,30], while at the macroscopic level one faces optimization problems different from the microscopic case, as shortly discussed in section 3.

The Sakaguchi model
We consider a system of N coupled oscillators, originally introduced by Sakaguchi [27], described by the Langevin equatioṅ where f i is an external constant force, and the Gaussian noise η i obeys the fluctuationdissipation relation Notice that we have chosen the system units such that the external force f i has dimension of frequency, which corresponds to taking the friction coefficient in eq. (2) equal to one. By introducing the complex order parameter where 0 ≤ σ(t) ≤ 1 measures the system coherence and ψ(t) is the common average phase, equation (1) becomeṡ where we understood the dependence on time. Let f 0 be the mean deterministic force, calculated over the N oscillator sample f 0 = j f j /N, we expect that the center of mass will oscillate with the frequency f 0 , so we can set ψ(t) = f 0 t + ψ 0 and thus we can redefine the dynamical variables as θ i = φ i − ψ(t), so as eq. (4) readṡ where we have redefined the external force as ω i = f i − f 0 . In principle eq. (5) represents a set of N coupled equations for the variables θ i , since σ is given by eq. (3). However, as N → ∞, one can replace the actual value of σ with its mean field value, and so eq. (5) becomes uncoupled. Such a mean field value can be obtained self consistently as discussed below. Eq. (5) corresponds to a Brownian particle moving in a periodic potential under the effect of a constant drift force ω i . Here and in the following we assume that the system reaches a steady state in the long time limit. In the course of this paper, we will discuss this assumption where relevant. The Langevin equation can be reformulated in terms of a Fokker-Planck (FP) equation for the probability distribution function (PDF) of finding the particle i at position θ at time t Thus, the stationary probability distribution function (PDF) of the position of such a particle reads [6,31] where I(x) = x 0 dy exp [β(Kσ cos y + ω i y)], and N is a normalization constant depending implicitly on β = 1/T , K · σ and ω i . The steady-state probability current thus reads J ss = N , and the particle steady-state velocity reads As N → ∞, we can adopt a continuous description, where the constant forces acting on the oscillators are distributed according to the probability distribution g(f ) with mean value f 0 . By introducing the shifted force distribution the self-consistent equation for the modulus σ of the complex order parameter, characterizing the degree of order or coherence in the configuration of the variables θ i is then given by [27] which can be decomposed into its real and imaginary part σ sin(ψ 0 ) = dω g 0 (ω) 2π 0 dθ p(θ, ω) sin θ .
By assuming that the force distribution g(f ) is symmetric around f 0 , and noticing that p(θ, −ω) = p(−θ, ω), the imaginary part on the right-hand side of eq.(10) vanishes, and so one is left with whose solution provides the mean field value for σ. As discussed in [27], for N → ∞ this model exhibits a critical coupling strength K c , such that for K > K c the systems exhibits a dynamical phase transition with synchronization σ > 0, while the system is incoherent for K < K c , and each particle described by the coordinate θ i oscillates with its proper frequency ω i . Thus, for K K c we expect σ to be positive but small, and we can expand eq. (13) in powers of ǫ = Kσ/T , obtaining while expanding eq. (8) the average velocity of the dynamical variable θ reads It is worth to note that the first two coefficients of the expansion (14) were recently rederived in ref. [33] by using a different approach: the author mapped the deterministic time evolution of the Kuramoto model order parameter into a stochastic process as given by eq. (1), and applied a fluctuation relation to the thermodynamic irreversible work done on such a stochastic system. Inspection of eq. (14) provides the critical coupling strength for which a nonvanishing solution to that equation appears The value of the order parameter σ as a function of K and T , for K > K c can be obtained by solving eq. (14), which gives where we have introduced the quantity ∆K = K − K c , and .
The figures contained in this manuscript and discussed in the following, are made by using the approximate expressions (15) and (17), and choosing the values of the parameters such that σ 0.5.
A few comments on the velocity v θ (ω) are now in order. Such a quantity represents the velocity of the dynamical variable θ and is thus the velocity deviation of a particle under the effect of a force ω with respect to the center of mass velocity f 0 . Inspection of eq. (5) suggests that for K < K c (thus for σ = 0) v θ (ω) = ω. On the other hand v θ (ω) goes to zero as K increases above K c : the higher K the higher are the barriers of the periodic force in eq. (5), while σ is also an increasing function of K.

Macroscopic power grids
Here we want to establish contact between the network of microscopic oscillators discussed in the previous section, and a corresponding macroscopic model used to describe macroscopic power grids. This analogy will be useful in the next section, where we will introduce the thermodynamic forces and powers characterizing the microscopic model. Given that in the macroscopic realm, one deals with alternating current (AC) networks, in terms of equations, this amounts to write a system of equations for the phase angle φ i for both the generators and the users [30] which correspond to an extended KM with inertia coefficients M i (representing, e.g., the large rotational inertia in turbine generators), viscous damping D i , power injection ω i (consumption if ω i < 0) and power flows along the lines a ij sin(φ i − φ j ) with coupling strength a ij . This model exhibits a synchronized phase that corresponds to a power grid that operates in a steady state with spatially uniform frequency.
In this kind of power network one may want to determine the optimal operation conditions to, e.g., avoid blackouts. For example in a static grid with a fixed number of producers and users, where the generators have a maximum power capability, and the users are characterized by a well known average consumption, one may want to optimize the grid in order to avoid that the consumers' load exceed the generators capability, leading to blackouts. This can be done by using optimization algorithms for graphs, see, e.g., [34].
Another, perhaps more interesting optimization problem considers a power grid as a dynamical system, where both energy producers and users can be dynamically connected or disconnected over a given time period. On the producers side, this is the case of renewable energy sources, such as wind and solar power, which are stochastic in nature and often uncontrollable, resulting in severe difficulties in the maintenance of the balance between load and generation [35]. Thus in the future opportunistic users can access the energy system according to the availability of system resources and differently from the "always-on" demand of traditional energy users, their consumption can exhibit peaks of activity. In this scenario, the challenge is to coordinate and manage dynamically interacting power grid participants.

Stochastic Thermodynamics of the microscopic model
In section 2 we have discussed the dynamical properties of the model system we will use in the present paper. We can now turn our attention to its thermodynamic properties, namely the input and delivered power, and the system efficiency as a global motor.
One can consider two possible scenarios as far as the forces applied on each single particle are concerned.
In the first case the forces acting on the motor system can be either positive or negative, thus resembling the macroscopic power grids of power plants and consumers considered, e.g., in [36,37,30,38]. Differently from those works, we consider here microscopic oscillators, in the over-damped regime, and with white noise acting on them. In this scenario, taking inspiration from the macroscopic realm, one may call users those oscillators with a negative force acting on them f i < 0, and producers those oscillators with a positive force f i > 0, and a single force distribution characterizes the system.
The second possible scenario resembles the case of molecular motors, where both a negative (f − i < 0) and a positive force (f + i > 0) are applied on the same particle i. This is the case in, e.g., biological molecular motors such as kinesin and myosin [39,40] where the energy extracted by ATP hydrolysis drives the motor forward (corresponding to f + i > 0) while the motor does work to carry a cargo, modelled by a negative load (corresponding to f − i < 0) In this case one deals with two different distributions of forces, g + (f + ) and g − (f − ). The homogeneous case, where the same positive f + and negative force f − where applied on all the motors, modelled as diffusing particle on a lattice with an exclusion rule, was studied, e.g., in [8,22,23].
In both scenarios, in order for the system to behave globally as a motor, and to perform work against the negative forces, we must require the center of mass to have an average positive velocity, and thus f 0 > 0.
In the following, we will consider, for both scenarios, the delivered power P out and the input power P in , and optimize P out wrt different parameters. We will thus characterize the efficiency at maximum power (EMP) η * = P * out /P * in mainly in proximity of the dynamical phase transition, and where possible, we will discuss the thermodynamic properties of the system in the whole range of parameter space.

Single force distribution
We consider here the stochastic thermodynamics of a system with either negative or positive forces applied on each oscillators, and distributed according to the single PDF g(f ).
We can thus introduce the relevant thermodynamic quantities, namely the average input power, absorbed by the producers, and the average output power released by the users. Recalling that v θ as given by eq. (8) gives the deviation of the i-th particle's average velocity from the center of mass velocity f 0 , the average output and input power read while the thermodynamic efficiency of the system reads Substituting eqs. (15) and (17) into (22) and (23) the output and input power becomes, up to the second order in K − K c , where with analogous definitions for P > 0 , I > 2 and I > 4 . We notice that, in absence of partial synchronization (K < K c , σ = 0), i.e., when the users and the producers are decoupled, eq. (22), (25), and (27) predict that the delivered power is negative. Since v θ (ω) = ω for K < K c , the users oscillates with their proper frequency (force) which is negative, and so the product of the applied forces times the average velocity is positive. The term v θ (ω) in eq. (22) is always negative, as the integration variable runs over negative value. Thus by increasing K above K c the modulus of v θ (ω) decreases and tends to zero for very large K. This implies that for some value of K the rhs of eq. (22) becomes positive, such values depending on f 0 and T , and on the details of the distribution g 0 (ω), e.g. its width.

Optimization
Here we aim at optimizing the delivered power eq. (22) wrt some of the relevant parameters.
Optimization wrt the coupling strength ∂P out /∂K = 0 gives: Recalling that the average velocity v θ goes to zero as K increases above K c , we find that, for ω < 0, v θ (ω, K) is an increasing function of K, ranging from ω for K < K c and approaching zero as K → ∞. Thus, from eq. (30) it follows Similarly one finds Thus, if K is the free parameter, the optimal delivered power is achieved for K → ∞, corresponding to the limit of strong coupling between users and producers, with an EMP η * = f − / f + as obtained by eqs. (22) and (23), where f − and f + are the average negative and positive forces, respectively.
No similar inequalities can be found when one tries to maximize P out with respect to other parameters, for example f 0 . So one should consider specific cases for the force distribution in order to study the relevant thermodynamic quantities.

A specific distribution
In order to exemplify the results discussed in this section, here we consider the specific distribution where s 2 is the variance of the distribution, with s > f 0 > 0, i.e., there are just two types of oscillator, the users with an applied force f 0 − s < 0 and the producers with an applied force f 0 + s > 0. The shifted force distribution thus reads For such a distribution the critical coupling strength reads This corresponds to the bimodal distribution considered in [41], where the linear stability of the incoherent solution p(θ, ω) = 1/2π of the FP equation (6) was studied, corresponding to the non-synchronized phase σ = 0. The authors showed that in the limit N → ∞ the FP equation eq. (6) exhibits a steady state solution for s < T and K > K c , while for s > T and K > 4T the system exhibits an oscillatory state. In the following we will thus take s < T , and use the steady state solution (7) to eq. (6).

Equations (25)-(26) thus become
We recall that for K < K c (i.e. for σ = 0), v θ (s) = s, and because of the condition s > f 0 we have P out < 0, i.e. when the producers and users are not coupled, the users oscillates with their proper frequency f 0 − s resulting in a negative P out . The delivered power will become positive for some value of K > K c , when v θ (s) in eq. (36) becomes smaller than f 0 .

Optimization
Here we optimize the delivered power for the force distribution (34), which corresponds to a system where on each oscillator there is either a positive f 0 + s or a negative f 0 − s force with probability 1/2. From eq. (15) we easily obtain the expression for the velocity deviation from the center of mass up to the fourth order in σ v θ (s) = s 1 − while the order parameter, as given by eq. (17), becomes up to the order 3/2 in ∆K.
We can now optimize P out , as given by eq. (36), wrt to different parameters: i) By optimizing wrt to K at fixed f 0 and s: ∂ K P out = 0, one obtains since v θ (s) is a decreasing function of K, as already discussed above for a general force distribution g 0 (ω). ii) By optimizing P out wrt the average force, at fixed K and s, ∂ f 0 P out = 0, one obtains and since the condition s ≥ v θ (s) > 0 holds for any K > 0, we have s > f * 0 (s, K) > s/2. We can thus calculate the delivered and the input power, and the efficiency at the maximum where we have used (38) and (39) 43) and (44) suggests that, for fixed s, when K increases above K c , the optimal output power (42) increases, the optimal input power (43) decreases, and this results in an increase of the EMP (44). This is a consequence of the fact that v θ (s) → 0 in the limit K → ∞, where η * = 1/3. Inspection of figure 1, as well as of eq. (45), suggests that a higher degree of quenched disorder, as parametrized by s, leads to a larger EMP close to the dynamical phase transition. However, the analysis of the behaviour of eqs. (42), (43) and (44) at fixed ∆K and varying s is not so straightforward. Graphical analysis of eqs. (42), (43) (not shown) indicates that both P * out and P * in increase with s, with P * out increasing faster. This graphical check can be done in the range of parameters where eqs. (38) and (39) holds, i.e. close to the critical point. However, it is worth to note that for a fixed ∆K, one finds P * out (s = 0) = 0. Furthermore inspection of eq. (5) also suggests that v θ (s) → s, as s → ∞, and being P * out a positive quantity, it must have at least one maximum for s ∈ [0, +∞[. On the other hand, P * in , as given by eq. (43) in an increasing function of s. Accordingly η * has at least one maximum for s ∈ [0, +∞[. iii) In order to optimize P out , as given by eq. (36), with respect to the quenched disorder standard deviation s, one has to solve the equation ∂ s P out = 0. This equation has no analytic solution s * , but it can be solved numerically, in order to find the EMP for different values of f 0 and K, as shown in fig. 2. Still one has the constraint s * > f 0 , in order for the force on the user to be negative.
One can also find an approximate expression for s * at the lowest order in ∆K and s/T as follows.
By solving ∂ s P out = 0 for K < K c , one easily finds that the maximum is given by s * 0 = f 0 , and thus P * out = 0. For K > K c , one can derive P out wrt s, and then expand the equations up to the second order in (K − K c ) and ∆s = s − s * 0 , solving for s, and plugging the value s * that maximize P out into the expression for η, one finds where the last expression gives η * to the lowest order in ∆K and f 0 /T . It is worth noting that the coefficients in the series expansions are the same as in eq. (45). Note that ∆K depends implicitly on s * through K c , as given by eq. (35), when one replaces s with s * . Inspection of eq. (36) suggests that the optimal quenched disorder standard deviation s * increases as f 0 increases, but this in turn leads to a smaller value of ∆K for fixed K, as K c (s = s * ) also increases. Thus increasing f 0 , and optimizing P out wrt s drives the system towards smaller values of synchronization σ, resulting in a smaller η * . This is in agreement with the results reported in fig. 2, showing that the EMP close to the dynamical phase transition, is enhanced by decreasing the applied average force f 0 .

Gaussian distribution
We now consider the following distribution for the force From eq. (16) we can calculate the critical coupling strength and from eqs. (22) and (23) we can calculate the output and the input power for K < K c which read with lim s→0 P out,0 = 0, lim No analytic result can be obtained for P out , and P in and thus for the EMP when K > K c for the Gaussian force distribution (47), at variance with what has been done in the previous section. However, one can resort to numerical calculations, to integrate numerically eqs. (22)-(23), find the optimal value of P out wrt to some parameter, and thus calculate the EMP. This procedure has been followed in order to calculate the EMP as obtained by maximizing P out wrt to the mean force f 0 . The results are reported in fig. 3: while for small ∆K the EMP is larger the smaller the variance, for sufficiently large ∆K we find the same tendency observed for the bimodal distribution, namely η * as a function of ∆K is larger, the broader is the force distribution, and thus the degree of quenched disorder.

Distribution of positive and negative forces on the same particle
In this section, we consider the case where on the same particle, whose dynamics is described by eq. (4), two forces are applied, f i,− < 0 and f i,+ > 0, distributed with two PDFs g − (f − ) and g + (f + ). In this framework the output and input power read respectively.

No disorder
We start our analysis by considering the trivial case where the same forces f + and f − are applied on all the particles, with f 0 = f + + f − > 0 and f − < 0. So, the total force distribution reads g(ω) = δ(ω − f 0 ), and the delivered power (53) becomes which is independent of the coupling K. Thus by maximizing P out wrt f − one finds that the EMP is always η * = 1/2, as in the linear regime case [2,6,7,42].

Bimodal negative force distribution
In order to increase the complexity of our system, we consider here the following distributions so as the total force distribution reads We introduce the following variables Given the expression for the average force, we obtain a condition on x: The delivered power eq. (53) thus becomes while the input power reads It turns out that with this choice of the force distributions the order parameter σ and thus v θ depends only on the negative forces, indeed we have and g 0 (ω) is the function appearing in the self consistent eq. (10). By deriving P out with respect to the disorder degree parameter y, we obtain ∂ y P out = −(v θ (y)+yv ′ θ (y)) < 0, where we have assumed y > 0 without loss of generality, i.e. P out decreases monotonically with y. Thus the optimal P out is trivially obtained for y = 0 which corresponds to the case with no disorder, with η * = f − /f + .
We now optimize P out eq. (62) wrt the variable x, which is equal to the mean negative force, and find and thus and finally the EMP reads We notice that while for K < K c one finds v θ (y) = y, and therefore in the uncoupled regime η * < 1/2, for K > K c the velocity v θ (y) is a decreasing function of K, thus for a fixed y the coupling reduces the spread around the mean velocity f 0 , and thus in the limit of large K one recovers the single particle EMP η * = 1/2. Similarly, P * out (66) is an increasing function of K, because of the decreasing behaviour of v θ (y). By expanding v θ in powers of K, and noticing that the order parameter σ is given by eq. (39) , with the substitution s → y, we obtain to the lowest order in ∆K and y. We notice that the quantity y 2 is equal to the variance of the negative forces' distribution, thus the introduction of disorder in the force distribution, on the one hand reduces the EMP for an uncoupled system (K < K c , first term on the rhs of eq. (68)), on the other hand it increases the slope of the EMP above the critical coupling. In fig. 4 the EMP η * is plotted as a function of K for different values of y.

General case: optimization
We now consider the case where the system exhibits a distribution of both positive and negative forces g − (f − ) and g + (f + ). The distribution of the total forces on each particle thus reads and where y ± = f ± −f ± . The distribution g 0 (ω) is symmetric around ω = 0 if both the distributions g ± (f ± ) are symmetric around their respective average valuef ± , a symmetry that we assume in the following.
We have thus and similarly where in the last equality we have used the above mentioned symmetry of g ± (f ± ). Thus, if we want to optimize P out wrt tof − , we obtain and We can thus evaluate the optimal delivered power below the critical coupling, and for very large coupling constant Close to the critical K up to the second order in σ, by using eq. (15), we find P * in =f 2 + 2 + y 2 + + dy + dy − g + (y + )g − (y − ) y + (y − + y + )K 2 σ 2 2 [T 2 + (y − + y + ) 2 ] ≃f 2 + 2 + y 2 to the lowest order in ∆K/T . We obtain finally the EMP . (78) inspection of this last equation suggests that the maximal slope of η * as a function of ∆K is obtained for y 2 + = 0. As far as the variance of the negative forces is concerned, we find a similar scenario as in the previous section: while y 2 − reduces the EMP for the uncoupled system, above the critical coupling the EMP increases faster the larger is y 2 − .

Conclusions
In the present paper we have investigated the thermodynamic properties of a model of microscopic oscillators, subject to thermodynamic forces. We considered the effect of the disorder on the delivered and injected power and on the EMP, and discussed the critical behavior of such quantities for different force distributions. We considered two forces distribution types, one that resembles the macroscopic power grids, and one that resembles a system of interacting autonomous motors.
For the first type of force distribution we find that, at fixed coupling strength, a larger degree of disorder leads to an increase in the EMP, at least close to the critical point.
For the second type of force distribution we find that while the disorder reduces both the optimal P out and the EMP below the critical coupling, above the critical point the EMP rate as a function of ∆K increases as the degree of disorder increases.
Thus, ideally the system with the optimal thermodynamic performances is characterized by a strong coupling (K → ∞) or absence of force disorder. However, in a real system one may have to deal with a finite coupling strength and an intrinsic non-vanishing degree of disorder in the force distribution. The results contained in this paper characterize the thermodynamic properties of such systems.
While we were able to calculate the expansion of the energy rates and of the EMP close to the critical point, we found that the EMP does not exhibit any universal behaviour, at variance with the single device case. On the contrary, the results contained in this paper depends strongly on the details of the force distribution. For example, in eq. (78), one recovers the EMP value in the linear regime (η * = 1/2) only when the disorder vanishes.
One of the limitations of the present model is that it exhibits an "all-to-all" coupling, while in a real system the interaction can depend on the distance between the network nodes and on the network topology, thus one should replace the interaction strength K in eq. (1) with an interaction matrix K ij . The thermodynamic properties of this extended model are certainly worth to investigate.
Furthermore, the values of the EMP reported in figs. 1, 2, 3 are quite small, thus the characterization of the response of the network injected and delivered power, or of its efficiency to a change in the network topology is certainly worthy of future investigation. For example, one may want to find the connectivity matrix between the different nodes that optimize the relevant thermodynamic quantities.