Mesoscopic virial equation for nonequilibrium statistical mechanics

We derive a class of mesoscopic virial equations governing energy partition between conjugate position and momentum variables of individual degrees of freedom. They are shown to apply to a wide range of nonequilibrium steady states with stochastic (Langevin) and deterministic (Nos\'e--Hoover) dynamics, and to extend to collective modes for models of heat-conducting lattices. A generalised macroscopic virial theorem ensues upon summation over all degrees of freedom. This theorem allows for the derivation of nonequilibrium state equations that involve dissipative heat flows on the same footing with state variables, as exemplified for inertial Brownian motion with solid friction and overdamped active Brownian particles subject to inhomogeneous pressure.


Introduction
From equilibrium statistical mechanics we are accustomed to the idea that there is energy equipartition among all quadratic degrees of freedom of classical systems, and that the "energy bit" corresponds to k B T /2, half of the temperature times the Boltzmann constant. While momenta usually appear with the quadratic contribution of the kinetic energy in the Hamiltonian H, for a position variable q i one has more generally that it is the average of q i ∂ q i H which equals the energy bit. The sum over all degrees of freedom yields the virial theorem [1,2], which connects the average total kinetic energy with the term i q i ∂ q i H named virial by Clausius.
Out of equilibrium, the equipartition of energy is not granted. Indeed, recent experiments with heat-conducting metals show intriguing deviations from equipartition, related to enhancements of low-frequency vibrational modes that may become even "hotter" than the highest boundary temperature [3].
Similar deviations from equipartition are observed for strongly heated cantilevers [4] and Brownian particles [5,6]. These are some out of many manifestations of nontrivial effects characterizing systems driven far from thermodynamic equilibrium. They imply the need for a critical revision of statistical mechanics equilibrium results, with the aim of pointing out possible extensions to nonequilibrium conditions.
In this work we discuss a generalization of the equipartition theorem, formulated in the context of modern nonequilibrium physics. It takes the form of mesoscopic virial equations (MVEs), involving kinetic and dynamical aspects specific to pairs of momentum-position conjugate variables. A MVE determines how thermal energy is distributed between any such pair of variables. For Langevin dynamics, we discuss both the inertial and the overdamped versions of the equation; the former is easily extended to cover Nosé-Hoover dynamics for thermostated simulations. Summation of a MVE over all degrees of freedom generates the virial theorem, which we discuss also for the case of explicitly nonconservative forces. That the virial theorem holds beyond thermal equilibrium should not come as a surprise, since it is a result derivable in classical mechanics without appealing to statistical arguments ‡. The simple mathematical derivations we employ are slightly different from the conventional line of arguments dating back to Chandrasekhar's work [8,2]. The main novelty of our approach is that we work consistently in the context of nonequilibrium systems, and that our derivations easily carry over to deterministic thermostats. Moreover, we characterise energy partition even for collective macroscopic variables, such as single normal modes, out of equilibrium. We further show that our results allow for the derivation of equations of state for nonequilibrium steady states. As an illustrative example, we provide a full derivation of the pressure equation for a well-known model of active matter [9].

Langevin dynamics
Consider N interacting particles evolving in d dimensions, with generalised coordinates {q i , p i }, i = 1, . . . , N d. Each degree of freedom has mass m i and the total energy is given by the Hamiltonian where U ({q i }) contains a confining potential energy that allows the system to reach a stationary state in the absence of external, time-dependent driving. In addition, nonconservative forces f i could also be present. Each particle is coupled to a Langevin thermostat with damping constant γ i , so that the general equations of motion reaḋ Here, the ξ i represent Gaussian white noise with correlation We first consider the case of independent heat baths in local equilibrium at temperature T i , for which the fluctuation-dissipation theorem implies a diagonal diffusivity matrix In section 6 we will show an example of a non-diagonal temperature matrix emerging for the normal modes of coupled oscillators. Note that a spacedependent noise is included in this formalism, since T i may be a continuous function T (q i ) of the particle position.
Our results directly derive from the formula for the time derivative of the average of any state observable O(t), where L is the backward generator of the dynamics. For the Langevin equations (2) it can be derived with Itô's formula [10] and is given by For example, a set of relations emerges from the position-momentum observable O = p i q i §. Plugging it into (3), we obtain d dt § This derivation is formally identical to the one employed in most quantum mechanics textbooks, e.g. [11].

Using then
and removing all time derivatives by the assumption of stationarity, this is turned into a MVE for the conjugated pairs q i , p i : The virial theorem follows by applying N d i=1 to both sides of (6). The equipartition theorem is recovered in equilibrium (f i = 0, T i = T ∀i), where averages may be performed with the Boltzmann weight exp(− H k B T ) and all terms in (6) are equal to k B T .
As a basic exemplification, consider a unit-mass particle moving in the q = (x, y) plane, subjected to the potential U (x, y) = x 2 − xy + y 2 and to the nonconservative shear force f = α y e x (e x being the x-axis versor) . In Fig. 1(a) we display a numerical validation of the MVE (6). Note that energy equipartition -with virial and (twice) kinetic contributions amounting to k B T -is achieved only in equilibrium (for α = 0). This example also illustrates that specific degrees of freedom can still be cooled down under nonequilibrium conditions (for α = 0) [12], despite energy being constantly supplied to the particle.

Nonequilibrium equations of state
Throughout the text we employ α as a dimensional constant measuring the departure from equilibrium.
m i is recognised as the average heat flow into the i-th reservoir, and in a steady state one gets the Harada-Sasa formula [13,14] Combining now the MVE (6) with (8), we find If the system is in thermal equilibrium, then J i = 0 ∀i, and (9) constitutes the standard starting point for deriving equations of state. For N interacting particles within a container of volume V it is useful to separate the contribution of the external conservative forces F ext (comprising confining wall forces F w , gravity, etc.) from that of the inter-particle interactions F int . The sum over all degrees of freedom of ∂ q i H q i gives both the external virial as a function of the pressure P , . Under equilibrium conditions, from (9) thus descends which can for example be used to derive the van der Waals equation [15]. Since the above deduction only relies on the homogeneity and isotropy of the system, it is directly applicable to nonequilibrium situations that enjoy these symmetries. We may think about systems with equal particles and homogeneous dissipation (T i = T and γ i = γ ∀i), say. As in the case of Fig. 1(a), nonequilibrium steady states are maintained by the action of the nonconservative forces, which contribute in (9) the additional nonequilibrium virial term C ne ≡ − N i=1 f i q i . Different cases should be distinguished, depending on the nature of f i . If f i is an external driving, such as the shear force of section 2, C ne combines with the conservative external forces to produce the pressure, Hence the equation of state (11) is generalised to where Q * = i J i /γ is the total heat flowing into the reservoir during a characteristic time 1/γ. On the other hand, if f i is a dissipative interaction force between particles (14) for a single particle (N = 1) at temperature T = 1 (natural dimensionless units). The negative sign of the heat flowing into the reservoir, Q * < 0, is consistent with a positive heat absorbed by the system when solid friction dissipates energy (α > 0).
(e.g. describing binary inelastic collisions in granular gases [16]), then (10) holds true and the nonequilibrium virial C ne figures explicitly in the equation of state Notice that these equations of state include not only equilibrium thermodynamic variables but also an unusual heat-flow contribution Q * = i f iqi /γ, which stems solely from the nonconservative driving because stationarity implies As a simple illustration of the role of Q * , consider N independent particles with unitary mass, again in the xy-plane. Each particle is subjected to a Langevin bath of uniform temperature T , to a confining potential U w (x, y) = 1 12 (x 12 + y 12 ) (thus F w = −∇U w ), and to an additional solid friction f = −α v/|v| of constant magnitude α ≥ 0 [17,18]. In the presence of this non-conservative friction, a steady state is generated in which heat is continuously taken from the Langevin bath and delivered to the substrate. However, the symmetry of the problem implies that C ne is zero. In view of the particles' mutual independency, also C int is exactly zero, and each of the remaining terms in (14) amounts to N times the single-particle contribution. In Fig. 2 we display each term in (14) for single-particle simulations of the system at various α. As depicted, Q * = 0 in equilibrium (α = 0), while Q * is negative in nonequilibrium (α > 0) and gives important contributions that guarantee the validity of the equation of state.

Overdamped dynamics
If one considers time scales much larger than the characteristic relaxation times of momenta, i.e. γ i dt → ∞ [10], then J i /γ i → 0 and (9) reduces to the overdamped This corresponds to (6) after the substitution p 2 i /m i → k B T i , as it should be expected, since momentum is instantaneously thermalised by its own thermal bath in the overdamped limit. Of course, this relation can be derived directly by taking the overdamped limit of the diffusion equations (2): where i is the appropriate observable to plug in (3) to retrieve (15).
These results hold under the assumption that the dissipative force f i acts effectively on time scales much longer than 1/γ i . If instead f i is of order O(γ i ), energy dissipation interferes with the thermalization process of momenta, so that and thus yields an overdamped MVE which features nonequilibrium corrections to the bath temperature, of the form Active Brownian particles (see more details in the next section) can be taken as another example. In the overdamped limit, they are often modeled as colloidal particles driven by a propulsion force f p,i that is counterbalanced by an associated viscous drag force −α i p i . Together they combine into the non-equilibirium force f i = −α i p i + f p,i . If the friction forces are comparable in magnitude, that is α i /γ i = const in the limit γ i → ∞, equation (8) in the overdamped limit reads which implies a renormalised temperature for the overdamped MVE ¶ To avoid the issues related to the interpretation of the overdamped stochastic equations hereafter we consider additive noise only.

Overdamped active matter
Active Brownian particles are often employed as an overdamped model for the collective behaviour of motile bacteria and self-propelled colloids [19]. Their phase behaviour is currently much studied [20,21,9,22,23]. In this regard, the utility of the virial theorem was pointed out in Ref. [24]. Here we fully exploit the generalised virial theorem and show how our approach leads to a pressure equation for active particles confined by hard walls of arbitrary geometry.
We describe an ensemble of identical active Brownian spheres moving in a twodimensional volume V in terms of their positions r i = (x i , y i ) and velocity orientations θ i (hence, {q} = {r, θ}). Their overdamped equations of motion arė The active velocity of modulus v 0 is directed along the unit vector u(θ i ) = (cos θ i , sin θ i ), and can be formally interpreted as another realization of the nonconservative force f that breaks detailed balance. Specifically, f = −αp + f p (see section 4) so that v 0 u(θ i ) = µf p with µ = [m(γ + α)] −1 . Each particle experiences the others through the two-body force F int . No special symmetry is assumed for the confining hard walls acting via F w (r i ) at the container surface S. The Gaussian translational noiseξ and the choice of the observable O = r 2 i in (3) yields the overdamped MVE Summing over all particles, the term containing F w becomes the external virial which is most conveniently expressed in terms of the local particle density, ρ(r) = N i δ(r − r i ) , as Since the local stress tensor P is defined by the steady-state equation expressing momentum conservation [25], an integration by parts of (24) HereP V is the volume-averaged pressure, defined through the trace of the pressure tensorP V ≡ 1 2V V drTrP (r). One expects the pressure to be non-uniform due to particle aggregation at the boundaries [26,27] and phase separation [28] in the presence of activity, unless highly symmetric geometries are considered [9]. Note that in the momentum balance (25) the only external force is the wall interaction. Consistently with the assumption of a constant active speed v 0 , the self-propulsion force and the corresponding fluid friction are taken to balance each other, hence are not included in the righthand side of (25).
For the special case of hard walls, the external virial is proportional to the surface- [29,30]. Moreover, interparticle interactions do not contribute to the momentum flux at the wall, so that the surface-averaged pressureP S can only have a kinetic contribution [30,31],P S = k B Tρ S . The latter is an equilibrium result that can be employed here, since in this overdamped description momenta are assumed to be thermalised at the temperature T . This follows from the choice of the translational noise's correlation. Therefore one arrives at the important result that the external virial gives the mean force per unit area exerted on the container, In the bulk, the interaction term in (23) gives a contribution analogous to the corrections to the ideal gas pressure in an equilibrium system. Indeed, for large N , where F int = −∇U int , r ≡ |r − r |, and g is the nonequilibrium pair density correlation function. In general, g cannot be reduced to a function of the relative pair position, since the system is inhomogeneous [32]. The explicit nonequilibrium contribution in (23) (the term containing v 0 ) gives rise to the so-called swim pressure [26,28]. Using (3), this time with O = r i · u(θ i ), and summing, we readily obtain v 0D The first average on the r.h.s. involves the particle polarization at the wall, while the second one represents the correlation between interactions and polarization. The constant term v 2 0 is an enhancement of the kinetic "ideal gas" contribution due to the particles' activity. Putting everything together, we obtain the equation of statē This result is valid irrespective of the confining geometry, thus extending the results of [24] and substantiating the numerical evidence for the equality of (average) wall and bulk pressure [26,27]. It is worth mentioning that it is currently under debate [23,22] whether (29) deserves the name of state equation, given that the active particle pressure appears to depend explicitly on the interactions with the boundaries and not only on bulk properties (temperature, density, etc.) as it would happen in equilibrium.

Normal modes of coupled oscillators
The derivation of the MVE does not rely on the diagonality of the the matrix D ij , that is (6) also holds for systems in which the noise components are cross-correlated.
An instance of such a situation is offered by the analysis of the normal modes of a system with local reservoirs. For harmonic lattices [7], depending on the details of the forcing and on boundary conditions, the energy stored in long wavelength vibrational modes may be either enhanced or reduced compared to the average. Here, we illustrate the MVE in modes' space for a one-dimensional chain of N point masses coupled with quadratic-quartic interactions, thus going beyond the harmonic approximation. The stochastic equation of the normal modes, obtained by applying a linear transformation to the equation (2) for the oscillators' position and velocity [33], is where ω 2 k is the squared eigenfrequency of the k-th mode and B klrs is a tensor that emerges from the quartic interactions. The noise terms η k , are mutually correlated, as quantified by a symmetric matrix T kl of mode temperatures [12], which is in general not diagonal unless the system is in equilibrium. Without the anharmonic coupling, = 0, the average kinetic and potential energy of the modes satisfy where the first equality is analogous to (6), and the second amounts to (9) specialised to the present analysis. Notice that the kinetic and potential energy coincide for a given mode, but they in general differ for different modes, thus breaking full equipartition. With = 0 the modes' dynamics is coupled via the tensor B klrs and the MVE (6) becomes containing no explicit sign of the non-diagonal T kl , as anticipated above. Similarly, the heat-flux equation (8) becomes This represents the perfect starting point to study perturbative corrections to mode energies, given the Gaussian statistics of the X k 's for = 0. The last term disappears in equilibrium (T i = T ∀i) where the modes' position and velocity are on average uncorrelated, so that the equipartition for velocities Ẋ 2 k = k B T results from (34). Differently, the non-zero heat flux present under non-equilibrium conditions modifies the mode kinetic energy in (34). For small we can expand (34) as Here we used the symmetry of the tensor B together with Wick's theorem to break up the Gaussian correlations . . . =0 evaluated in the harmonic system [12]. An illustration of (35) is provided in Fig. 3(a) for a one-dimensional lattice with fixed boundaries immersed in a linear temperature profile. For = 0, due to the symmetry in the T i 's and in the boundary conditions, the modes enjoy a peculiar full energy equipartition [12] at the average temperature The anharmonic corrections allow energy to leak into the higher, more localised modes. The same qualitative behaviour is found numerically for increasing values of (Fig. 3b). This suggests that the energy repartition among modes is robust against the introduction of non-linearities and fairly well approximated by a first order perturbative calculation. Note that the total energy is insensitive to , namely k=0 T kk ∀ , since the total heat flux appearing in (34) is identically zero under stationary conditions, thanks to the potential nature of the interaction: k,l,r,s

Deterministic thermostats
The relations derived above for stochastic inertial systems remain valid in the zeronoise limit, where the dynamics becomes deterministic. Stationarity is then ensured by coupling the system to suitably defined thermostats. Examples are Nosé-Hoover thermostats, where extra degrees of freedom act as frictional couplings for the physical ones [34]. Similarly to Langevin dynamics, they guarantee canonical thermalization in cases of uniform temperatures, and they lead to non-zero heat fluxes if different temperatures are imposed on different degrees of freedom of the system. For lattices of oscillators interacting only via conservative forces and coupled to Nosé-Hoover thermostats at various temperatures, the existence of local energy equipartition is a common assumption needed for the local definition of temperature [35]. So far, it has only been observed in simulations for the masses not directly driven by Nosé-Hoover thermostats [36]. Here we provide a formal proof. We consider statistical averages with  respect to the invariant density, which, in general, may or may not coincide with time averages. Equality is in fact assured by the use of Nosé-Hoover chains of thermostats [37].
The Nosé-Hoover dynamics for unit masses is given bẏ where Θ i is an indicator function, which is 1 or 0 depending on whether the mass i is coupled or not to a thermostat. The auxiliary feedback variable ζ i aims at thermalizing p i at the temperature T i on a timescale τ . The backward generator associated to (37) is Following the scheme outlined above, we find the generalised MVE which includes the formal justification for the mentioned numerical observation of local energy equipartition if restricted to masses without a local thermostat [36], corresponding to Θ i = 0. The term stemming from the thermostat's force (that can be seen as another realization of the non-conservative force f i ), is identically zero only in equilibrium, where momentum and position are uncorrelated and p 2 i = k B T i holds also for the degrees of freedom coupled to thermostats.

Conclusions
For a wide class of nonequilibrium systems in steady states, including stochastic and deterministic thermostated dynamics, we have shown that the kinetic energy of a given degree of freedom is on average equal to the corresponding virial of the forces. An integration over all degrees of freedom of such MVE yields the standard (macroscopic) virial theorem and a variety of useful results for general nonequilibrium systems. It is indeed possible to follow the path valid for equilibrium systems, using the virial theorem as a tool for the derivation of equations of state that involve pressure, temperature and other observables. For inertial systems with dissipative dynamics, this leads to an intriguing relation between the virial, the temperature of the heat baths, and the heat flux into them. Similarly, for active Brownian particles the virial theorem represents a powerful tool for deducing an equation of state valid for arbitrary container geometries. A direct experimental verification of the fundamental mesoscopic virial relations (underlying all these results) would therefore be desirable. In boundary driven systems with conservative internal forces, such verification amounts to checking energy equipartition between momentum-position type conjugate variables.