Dissipation Function: Nonequilibrium Physics and Dynamical Systems

An exact response theory has recently been developed within the field of Nonequilibrium Molecular Dynamics. Its main ingredient is known as the Dissipation Function, Ω. This quantity determines nonequilbrium properties like thermodynamic potentials do with equilibrium states. In particular, Ω can be used to determine the exact response of particle systems obeying classical mechanical laws, subjected to perturbations of arbitrary size. Under certain conditions, it can also be used to express the response of a single system, in contrast to the standard response theory, which concerns ensembles of identical systems. The dimensions of Ω are those of a rate, hence Ω can be associated with the entropy production rate, provided local thermodynamic equilibrium holds. When this is not the case for a particle system, or generic dynamical systems are considered, Ω can equally be defined, and it yields formal, thermodynamic-like, relations. While such relations may have no physical content, they may still constitute interesting characterizations of the relevant dynamics. Moreover, such a formal approach turns physically relevant, because it allows a deeper analysis of Ω and of response theory than possible in case of fully fledged physical models. Here, we investigate the relation between linear and exact response, pointing out conditions for the validity of the response theory, as well as difficulties and opportunities for the physical interpretation of certain formal results.


Introduction
Since the 1980s, a dynamical systems approach to nonequilibrium phenomena has been developed, side by side with the increase of computing power and the ability to perform Molecular Dynamics (MD) simulations. Classical MD consists of algorithms meant to solve Newton's dynamical equations for interacting particle systems, that may be subjected to external fields. Such investigations result quite informative and practically useful, when quantum effects are irrelevant or negligible, which is the case of fluids under standard conditions.
Provided no dissipative forces are present, one speaks of equilibrium MD, which is widely used together with statistical mechanical relations, in order to compute: rheological properties, polymers behaviours in porous media, defects in crystals, friction between surfaces, formation of atomic clusters, features of biological macromolecules, formation of drugs, etc. The results are most satisfactory, and, in practice, difficulties only arise if the interatomic forces are too complicated to As a matter of fact, NEMD has become popular for its effectiveness in computing transport coefficients, such as the viscosity of a fluid [6]. It was also realized that such coefficients can be obtained from dynamical systems properties such as the sum of Lyapunov exponents or, under special conditions, just from the sum of the largest and of the smallest exponent [6]. In turn, −Λ, the opposite of the phase space volumes variation rate Λ, which turns out to be proportional to the thermostatting (fictitious or synthetic) term may also be proportional to the entropy production σ when thermodynamics applies, see e.g., Ref. [19]. (A similar idea had previously been independently investigated by Andrej, who followed a different approach [20]). This idea became very popular when Evans, Cohen and Morriss applied a representation of SRB measures to the fluctuations of Λ in a Gaussian isoenergetic shearing system [21], as in that case such a quantity equals the entropy production, when local thermodynamic equilibrium (LTE) holds. That paper derived and tested the Fluctuation Relation (FR), one of the rare exact relations for nonequilibrium thermodynamics based on microscopic deterministic dynamics. For a time reversal invariant model of shearing fluids, the FR states that positive values of the energy dissipation, are exponentially more probable than their opposite. This was interpreted as an explanation of the second law of thermodynamics, and motivated a flurry of activity. Moreover, it led various authors to identify the entropy production with the phase space contraction [22,23].
Although Λ is relatively easy to handle in dynamical systems theory, it was pointed out that its identification with σ is only justified in a limited set of cases [24,25]. Furthermore, in time dependent settings and in presence of fluctuations, the connection of Λ with physically measurable properties results rather loose. It was then shown that the physically relevant quantity in the cases described by NEMD is the so called Dissipation Function Ω, which preserves the meaning of energy dissipation rate even when LTE is violated, and thermodynamic quantities such as σ do not exist. Only in a few cases Ω equals −Λ, although steady state averages of Ω and −Λ may be simply related to each other.
Later Ω became the basis for a general response theory, in which it plays a role as fundamental for nonequilibrium states as thermodynamic potentials do for equilibrium states. It can be used to solve in terms of physically relevant quantities the Liouville equation and, consequently, to investigate the exact average response of the observables of ensembles of systems. Furthermore, Ω can be used to provide conditions for the evolution of single systems, [26,27]. (This should be contrasted with the fact that almost invariably response of systems to perturbations is given in terms of ensembles, and that a comparatively limited amount of works have tackled response beyond its linear approximations). Most importantly, these results are expressed in terms of a physically measurable quantity, rather than in terms of an abstract phase space quantity.
Unfortunately, it is rather complicated to thoroughly investigate the properties of the dissipation function in realistic particle models. Therefore, we consider mathematically amenable dynamical systems, even if their physical relevance may be limited, as often and successfully done in statistical physics, with the purpose of illustrating aspects of the relation between dynamics and thermodynamics; think for example of Helmoltz's celebrated theorem [28]. Recently, Qian, Wang and Yi have also followed this route, scrutinizing a number of interesting examples and investigating them in terms of the phase space contraction rate −Λ [29]. The result is a set of intriguing relations between different dynamical quantities and formal thermodynamic quantities. In consideration of the fact that small systems, hence non-thermodynamic fluctuating quantities, as well as time dependent situations are ever more relevant in science and technology, we are going to perform a similar analysis within the dissipation function formalism. This paper is organized as follows. Section 2 summarizes the theory of the Dissipation Function, also analyzing exactly solvable examples and numerical simulations. Section 3 compares linear Green-Kubo and exact response. Section 4 follows in part Ref. [29], using Ω rather than Λ. The two approaches appear as largely equivalent, since they both arise from the solution of the Liouville equation. Section 5 explicitly solves the cases of simple attractors, and of small dissipation, in the sense of small Λ. Section 6 discusses with examples the conditions under which the exact response theory based on Ω holds, showing that Ω has to be continuous at least, and that the initial phase space probability density f 0 has to be positive in all the space that can be explored by the dynamics. This condition, called ergodic consistency [30] is reminiscent, but weaker than the requirement that initial and asymptotic measure be absolutely continuous with respect to each other, also called equivalence [31]. If that is not the case, modifications to the response formulae are required [32].
Finally, Section 7 summarises our work, highlighting two points in particular: investigating simple systems allows us to understand properties of Ω that are not immediately evident in realistic particle systems; on the other hand we propose a kind of thermodynamic formalism for the analysis of dynamical systems, cf. Ref. [33], that is based on the language of nonequilibrium, rather than of equilibrium, statistical mechanics. Section 7 also stresses the role of Ω as the energy dissipation, in the case of interacting particle systems.
Remark: one of the anonymous referees pointed out that the scope of our investigation can be substantially extended, because: "most quantum transport effects could be recast into a classical simulation picture of nonlocal collisions", cf. [34,35].

The Dissipation Function
Consider a system whose "microscopic" state is described by a dynamical systeṁ Let (1) generate a flow on M, i.e., a map S t : M → M, such that S t (x), t ∈ R is the solution of (1) at time t with initial condition x. Further properties of M and V(x) will be specified when needed. Let µ 0 be a probability measure absolutely continuous w.r.t. the Lebesgue measure dx, and denote by f 0 (x) its density function, i.e., dµ 0 ( for any measurable set E ⊂ M and for all t ∈ R. If µ 0 is not invariant, it evolves under the dynamics, so that at time t it is expressed by µ t (E) = µ 0 (S −t E). Provided the flow is sufficiently smooth, µ t is absolutely continuous, and has density f t (x) that satisfies [27]: where Ω f t (x) is the dissipation function, defined by: and denotes the phase space variation rate. Its negative, −Λ(x), is called phase space contraction rate. In the case the dynamical systems represents an equilibrium system, with equilibrium distribution µ 0 , so that f t (x) ≡ f 0 (x), the dissipation function identically vanishes, Ω f 0 (x) ≡ 0. It is obviously true also the converse: if Ω f 0 (x) ≡ 0 then f 0 (x) is an invariant density. In the general case, the solution of (2), is given by [27]: where, for every phase function O we use the notation Given any observable O : M → R, we define its ensemble average according to the measure Under suitable conditions, discussed in Section 6, this average can be computed using the dissipation function as follows [27]: Equation (8) gives the response of an equilibrium system, characterized by the distribution µ 0 , invariant for the dynamicsẋ = V 0 (x), which at time 0 is perturbed so that its dynamics becomes Equation (1). The measure µ 0 is not invariant under the new dynamics, hence it evolves together with the averages representing the observables. At time t the initial measure has turned into µ t , and the initial average O 0 into O t .
In the long time limit, the response may be given by a non equilibrium invariant distribution. As (8) is exact, the existence of such an asymptotic average is equivalent, as a necessary and sufficient condition called ΩT-mixing, to the convergence of the following limit: for any observable O. A sufficient condition for the existence of L O is [26]: A further useful asymptotic property is obtained by considering two observables O and Q. In this case, we say that the system is T-mixing [26] if We observe that setting Q = Ω f 0 in (11) we obtain a weak form of ΩT-mixing, namely given that one has, cf. Section 6 and [27]: Obviously, T-mixing implies weak ΩT-mixing. We observe that these properties may hold on a specific set of observables for a given f 0 but not for other initial distributions.
A fundamental difference between T-mixing and the mixing property of standard ergodic theory is that T-mixing concerns a the decay of correlations w.r.t. the initial state µ 0 , hence it may properly describe the irreversible relaxation of a given system to a stationary state. Mixing, instead, corresponds to a decay of correlations within a given steady state, since it refers to an invariant measure [5]. We refer to Refs. [5,27] for a broader discussion of these issues.

Examples
The physical meaning of the Dissipation Function can be appreciated by considering some special choices of vector field V, and of initial probability density f 0 . First, consider the NEMD systems for which Ω has been introduced, to see that in such cases it is directly related to the instantaneous dissipative flux. The simplest of these systems have the following equations of motion: where Γ = (q, p) collectively represents all coordinates q and momenta p of the N particles, and F e is the external driving field, coupled to the particles via the "charges" C i (Γ) and D i (Γ). The term α is a Lagrange multiplier meant to implement some constraint, so that the system reaches a steady state despite the driving field. In the NEMD literature α is said to implement a deterministic, time reversible thermostat [6]. If this thermostat is set to 0, one obtains the adiabatic equations of motion: which imply a continuous growth of energy. Then, one must assume that the support of f 0 is densely explored by a single equilibrium trajectory [30]. Denoting by Φ the internal interaction potential energy of the N particles, the dissipative flux, J, is obtained from the time derivative of the internal energy under the dynamics (15), which represents how much energy has to be dissipated in order to keep the steady state. The result is: where V is the volume of the system. Consider a few common examples.

Microcanonical Ensemble
For the isoenergetic dynamics, in which the internal energy is conserved, i.e., H 0 is fixed anḋ H 0 = 0, the equilibrium ensemble f 0 is microcanonical. Therefore, for an equilibrium Hamiltonian N particle system in d dimensions, one has [30]: where K is the kinetic energy of the N particles, while as in the original paper [36]. Then, one may write: As even the O(1/N) correction in Equation (18) is proportional to α, one observes that in this case both Ω and Λ are proportional to the rate of energy dissipation (17).

Gaussian Isokinetic System
In this case the kinetic energy K is fixed, and one may formally write K = (dN − d − 1)k B T = (dN − d − 1)/β, if also the center of mass is fixed, to prevent contributions from drift, so that: In turn, Λ is still expressed by (18), with a different O(1/N) term. In this case, f 0 writes: where Q is the normalization constant. The result for the dissipation function is: where the bar over the observables denotes time average. Combining this with Equations (18) and (19), one obtains where σ is the entropy production, if local thermodynamic equilibrium is verified. We observe that −Λ and Ω differ by a term proportional to the fluctuatingḢ. If these fluctuations average to zero, the averages of −Λ and Ω coincide, however, depending on how the respective fluctuations are distributed, various relations may apply to one and not to the other. For instance, the steady state fluctuation relation may hold for Ω but not for Λ, [37,38].

Nosé-Hoover Thermostat
In this case, there is one extra equation to be added to the equations of motion (14): where τ sets the relaxation time scale of the thermostat, and T is the imposed average temperature. Then, we may write: and for unperturbed Hamiltonian dynamics, one may write: In this case, f 0 is given by where M is the normalization. It follows that: which implies: We conclude that for equilibrium Hamiltonian dynamics, the driven nonequilibrium evolution is characterized by the dissipation function, which equals the entropy production, in case local thermodynamic equilibrium holds. Additionally, in this case −Λ and Ω differ by a term proportional toḢ 0 , and the same remarks concluding the previous subsection equally apply.

Canonical Ensemble for Generic Dynamics
The canonical ensemble describes the equilibrium state of a Hamiltonian particle systems in contact with a heat bath at temperature k B T = β −1 , with which it can exchange energy but not matter. Here, we consider a general case in which the vector field V(x) is not necessarily Hamiltonian, but there still is an energy function H. In addition, we assume that the initial distribution has the same form of the canonical Gibbs density: Z β being the normalization constant. Then, the Dissipation Function reads: where we use the notationḢ(x) = V(x) · ∇H(x). The Dissipation Function Ω f 0 is not necessarily zero, in this case, and the density evolves as follows, cf. (5): In the case of Hamiltonian dynamics, x = (q, p) ∈ R 2n , V(x) = J∇H(x) where J is a symplectic matrix, we haveḢ(x) = V(x) · ∇H(x) ≡ 0 and Λ(x) ≡ 0. Thus, in the equilibrium canonical ensemble, where (26) is invariant and the phase space volume is conserved, the dissipation function is identically zero The extended canonical equilibrium distribution of the Nosé-Hoover thermostatted system can be framed in this way: its dynamics includes a degree of freedom ζ that represents the heat bath, and system plus bath constitute a canonical equilibrium state. For instance, consider a system made of N identical noninteracting one-dimensional oscillators in contact with a bath at temperature T, implemented by a single Nosé-Hoover thermostat. This thermostat introduces an indirect coupling among the oscillators, which however does not represent a proper interaction leading to LTE. We illustrate these facts considering both linear and nonlinear oscillators, i.e., potential energies U(q) = −ω 2 0 q α+1 /(α + 1), with α = 1 and α = 3. The equations of motion are where T is the target kinetic energy per particle, is the kinetic energy per particle, and τ the relaxation time.
In order to elucidate the non-thermodynamic nature of the model (28), we simulated a system of N = 30 oscillators in the linear case, α = 1, with m = 1, ω 2 0 = 2, τ = 2 and T = 100. As expected because of the improper particles coupling expressed by ζ, we observe no equipartition; cf. Figure 1, which shows that the time average K i of the kinetic energies of single particles substantially change with their label i, while the time average K of the total is equal to T within numerical errors. The same happens also with α = 3.
To promote LTE in this system, we introduce nearest neighbour interactions as: in place of the second of (28), taking k > 0 as the coupling strength, and q 0 = q n , q n+1 = q 0 as boundary conditions. Unlike the previous non-interacting model, the linear (α = 1) and nonlinear (α = 3) cases behave in a different way. While still absent in the linear case, equipartition is established with the cubic local field. Figure 2 shows the time averages 1 J ∑ J j=1 p 2 i (t j ) of the kinetic energies of some particles sampled form a long trajectory of the coupled system with N = 10. We consider the cases with k = 1, 2, 5 in order to make evident the role of the coupling. These simulations show that increasing the particles coupling k speeds up the equipartition process. Changing the oscillators frequency, instead, has no effect on equipartition, cf. Figure 1. Now, focusing on the energy i of the system of oscillators alone, i.e., without the bath, and taking f 0 as in order to define a kind of reduced dissipation function for the dynamics with the bath, one obtains: which is not identically zero even in the equilibrium state. There is no immediate physical interpretation of the quantity Ω f 0 0 , but it clearly shows that, in this framework, the thermostat has to be included to properly speak of equilibrium.
, , , , of the kinetic energies of three particles (i = 1, 4, 8) sampled at t j = 10 3 · j, for increasing J. In the left panel, where J = 500, the outer continuous (red online) lines correspond to the coupling parameter k = 1, while the inner dotted lines (blue online) correspond to k = 2. In the right panel the focus is on a shorter interval, J = 150, where the relaxation phase can be better appreciated. In this case the outer lines (red online) correspond to k = 1, the middle lines (blue online) correspond to k = 2, the inner line (black) correspond to k = 5. Larger coupling enhances the convergence rate towards equipartition, as expected.

Hamiltonian and Gradient Vector Field with Canonical and Gaussian f 0
In the case of Hamiltonian dynamics perturbed by a gradient term [39]: where the unperturbed equilibrium density should be chosen consistently with the constraints that remain in place when the perturbation is switched on, in order to implement the above-mentioned ergodic consistency condition. For instance, if the system temperature remains fixed, one may take the canonical distribution, f 0 (x) = e −βH(x) /Z, that implies: One observes that Ω f 0 only depends on L if ∇L is orthogonal to ∇H. If ∇L is in addition a Hamiltonian perturbation, the perturbed steady state remains equal to the unperturbed equilibrium. Now, taking an initial uncorrelated Gaussian distribution: one obtains: J ij being the entries of J. Again, a Hamiltonian perturbation such that implies that the Gaussian (32) is invariant under the dynamics. These two examples constitute elementary warnings against too quick thermodynamic interpretations of generic particle systems, much as it may be interesting in certain circumstances [40].

The Dissipation Function and the Phase Space Contraction Rate
The Dissipation Function satisfies the Non Equilibrium Partition identity [27]: Using Jensen inequality, we can also write: obtaining: in accordance with both our physical intuition on the dissipation of energy in macroscopic systems and with the interpretation of this quantity as a Kullback-Leibler divergence between the initial and the evolved probability distribution in phase space [27]. One motivation to derive relations like these, see e.g., Refs. [27,30] where various others are derived, resides in the fact that they are only based on the dynamics and on the initial probability distribution, whatever phenomenon they represent, hence they hold very generally. At the same time, this also means that blind thermodynamic interpretations are to be avoided, although no difficulty arises from a purely dynamical systems standpoint.
The distinction between the dissipation function and the phase space contraction rate −Λ becomes particularly evident when considering particle systems. In fact, the steady state phase space average of −Λ equals the steady state entropy production of certain NEMD models [6,19], and it can be expressed as: µ being the corresponding invariant measure, and e p the steady state entropy production or the energy dissipation. On these grounds, Gallavotti and Cohen proposed −Λ as the steady state definition of entropy production tout court, [41], followed by Ruelle [23]. Earlier, a similar idea had been contemplated by Andrey [20], who referred to the variation in time of the Gibbs entropy. The Gibbs entropy, however, is not directly related, in general, to the entropy of a nonequilbrium system [42][43][44].
(Note: this statement does not concern the formally analogous definition of entropy common in the theory of stochastic processes, since they are not defined in phase space but in the space of observables.) Therefore, −Λ cannot be directly related, in general, to the entropy production, when considering time dependent states or steady state fluctuations, cf. Section 2.1 and Refs. [27,30,38]. On the other hand, it is interesting to investigate both −Λ and Ω as properties of dynamical systems, in a kind of extension of the thermodynamic formalism [33], that refers the non-equilibrium language.
Given an initial absolutely continuous probability measure dµ 0 (x) = f 0 (x)dx on the phase space M, and letting µ t be the evolution of µ 0 at time t, the corresponding average of the dissipation function is: which in the case of a uniform density f 0 yields: In this case, even Ω f 0 does not need to represent any entropy production rate: the equilibrium dynamics would have at least to be phase space volumes preserving, for Ω f 0 to be related to the dissipation of energy [30]. Using Equation (8) and the fact that Ω f 0 0 = 0, we may also write: In the case of a generic probability density f 0 , the exact response formula (8) and (40) yield: (41) which quantifies the difference between the mean phase space contraction rate and the energy dissipation in time. In the case ΩT-mixing, Equation (9), holds, t can be replaced by ∞ in (41), obtaining the asymptotic difference between the two quantities: From the definition of dissipation function Ω f 0 , i.e., Equation (3) with t = 0, one also observes that Therefore, the time integrated difference of −Λ and Ω f 0 merely amounts to a boundary term: which, however, may bear noticeably affect the statistic of fluctuations: the distribution of such boundary terms may, for instance, imply the validity of the FR for Ω f 0 and not for −Λ, cf. [37,38,45].

Exact Response and GK Formalism
In this section, we compare the classical linear response theory, or Green-Kubo formalism, with the exact response represented by Equation (8). Consider a perturbed vector field where may or may not be small. Then, the Dissipation Function can be written as: and the evolution of a probability density f 0 (x) is given by where S t denotes the flow of the vector field G (x). Thus, the exact response equation (8) reads In case G 0 derives from a Hamiltonian H 0 , the unperturbed dynamics preserves phase space volumes, and one can write: where L 0 = [H 0 , ·] is the Liouvillean corresponding to H 0 and [·, ·] denotes the associated Poisson brackets. Then, we may define the perturbation operator L p as: so that the (total) dissipation function can be written as: Substituting (51) in (48) and in case the initial distribution f 0 (x) is the equilibrium state of the unperturbed system, so that L 0 f 0 ≡ 0, one obtains: Provided the response Equation (8) applies, the formulae (52) and (53) for the evolution of phase space averages are both exact: there is no assumption on the size of and no approximation has been made. An approximate formula can be obtained in the small-limit, starting from Equation (47), by truncating to first order the exponential expansion, and recalling that Ω This yields: where the tilde symbol stresses that this is but an approximation of the probability density at time t, which is not guaranteed to be a probability density itself. Nevertheless, multiplying (54) by O and integrating over the phase space, one gets: In order to keep the truncation error bounded in time, one may assume that the time integral , which is reminiscent of the ΩT-mixing condition. The difference between the exact response formula, (48), and the approximate one, (55), amounts to: that, in general, does not vanish, and may grow large in time.
May the difference of responses Equation (56) become negligible, in some limit? Applying the transformation S −τ x = y to the second integrand of (56), one obtains: (57) where J S τ (y) = 1 since we are dealing with Hamiltonian dynamics-being J S τ (y) the Jacobian determinant. Thus, Equation (56) can be rewritten as: For every finite t, this quantity vanishes in the → 0 limit but, unlike the case of the Green-Kubo theory of linear response, the flow used in (58) is the perturbed one, S τ , not the unperturbed S τ 0 . As a matter of fact, the success of the Green-Kubo theory is based on its ability to compute the response by averaging with respect to the equilibrium ensemble for the unperturbed dynamics, and by simulating such an unperturbed dynamics [46]. The procedure, illustrated in many textbooks, including [46,47], uses the canonical ensemble for the unperturbed Hamiltonian H 0 , Z −1 0 exp{−βH 0 }, and assumes that the parameter is the magnitude of and external force F e . Within the framework of NEMD and Equation (14), one starts from the unthermostatted, α = 0 situation. Without loss of generality, the external driving F e can be written as a scalar quantity, since its direction can be embedded in the coefficients C i and D i . Then, the unthermostatted perturbed equations are derived from a Hamiltonian, H = H 0 + H p . Writing x = (q, p) to distinguish configuration from momenta, G p is obtained deriving H p : and the following holds: In this case, one has: which equals Equation (16), and shows again that Ω f 0 G , not −Λ, represents that energy dissipation rate. Now, following the Green-Kubo approach [46,47], which rests on the fact that one may use the unperturbed dynamics, and eventually write: to first order in = F e , for a generic observable O. The function R(t) = J O • S t 0 0 is called the response function. In particular, the difference in (56) between exact and linear response can be expressed as: That the full dynamics S τ can be replaced by the equilibrium dynamics S τ 0 is a far from trivial fact. In fact, van Kampen raised one of the best known objections to this theory, based on the observation that particles dynamics are typically chaotic [48]. He maintained that while there is no doubt that linear response works in practice, the derivations a la Green and Kubo do not explain why. He stated that: "Linear response theory does provide expressions for the phenomenological coefficients, but I assert that it arrives at these expressions by a mathematical exercise, rather than by describing the actual mechanism which is responsible for the response" [48]. It fact, to justify the Green-Kubo calculation, one must invoke non trivial ingredients, such as the large value of N, which makes the perturbation per degree of freedom negligible, and the mingling of contrasting contributions operated by chaos, see e.g., [45,49,50]. As such conditions actually hold in macroscopic particles systems under quite variegated conditions, the Green-Kubo reasoning explains the success of linear response.

Exact Response and GK Formalism: Example 1
Consider a particle in a two-dimensional space, whose unperturbed Hamiltonian is H 0 = (p 2 x + p 2 y )/2. Take the canonical ensemble f 0 = Z −1 0 exp(−βH 0 ) as the initial distribution. Introduce an external field with Hamiltonian H p = −Eq x , as a perturbation to obtain H = H 0 + H p . The corresponding equations of motion are: with the dissipation function Ω in which the momentum in x-direction plays the role of internal dissipative flux with E the corresponding thermodynamic force, i.e., J = −p x and F e = E. Note that the present example is an instance of the adiabatic system (15), in which C 1 = C 2 = 0 and D 1 = 1, D 2 = 0, see (59). Sticking to the notation of Section 2.1 in the following we denote by Γ = (q x , q y , p x , p y ) the phase space points and observe that the solution of (66) with initial condition Γ gives p x (Γ, t) = p x + Et and p y (Γ, t) = p y . Then, thanks to Equation (5), we compute the time evolution of the probability distribution function f t (Γ) starting from the integrated dissipation function: Following the discussion in Section 3, we expand the exponential in the limit of small external field E, and we neglect the terms of order higher than linear in E, to obtain the linear response for the evolving density:f The difference between exact and linear response can be further clarified by looking at the Liouville equation for f t (Γ) and proceeding à la Green-Kubo. Let us start introducing the Liouville operator of the perturbed system: where L 0 = p x ∂ q x + p y ∂ q y corresponds to the unperturbed system and ∆L ≡ L p = E∂ p x corresponds to the perturbation. Splitting the evolving ensemble distribution function as: the Liouville equation takes the form: Recalling that the equilibrium distribution f 0 (Γ) is invariant under the unperturbed dynamics-so that (∂ t + L 0 ) f 0 (Γ) = 0-we obtain the evolution equation for the deviation ∆ f t (Γ): Then, assuming that ∆L∆ f t (Γ) is negligible w.r.t. the other terms, we arrive at: Since (67) can be rewritten as where β f 0 = −∂ f 0 /∂H 0 , E = F e and p x = −J, according to standard notation. The solution to the evolution equation for the variation of the ensemble distribution function, is given by: where we have used the fact that the external perturbation is not time dependent.
Let us now compute the phase space average of an observable O, using f t (Γ) = f 0 (Γ) + ∆ f t (Γ). We can write: noting that S t 0 is the flow of the unperturbed dynamics. As an example, take O = p x , and compute its time evolving linear response. Since p x 0 = 0 and p 2 x 0 = β −1 , we have: Now, denoting by S t E the flow of the perturbed system, and recalling that p x • S t E = p x + Et and Ω f 0 E = βEp x , the exact response formula (8) yields: For such an observable, linear and exact response coincide. In general, we have: For an analytic observable that do not depend on q x and q y , we have: assuming that O = O(p x , p y ) is an analytic function, and denoting by O (k) the k−th partial derivative of O with respect to p x , at fixed p y . The difference between the exact and the linear response would read: Depending on the particular functional form of O this quantity may grow rapidly or less rapidly with t and, as long as O(S t E Γ) exists at all t, this difference at fixed t goes to 0 as E → 0.

Exact Response and GK Formalism: Example 2
Let the unperturbed Hamiltonian correspond to a harmonic oscillator H 0 = (p 2 + ω 2 q 2 )/2, and let the perturbation take the form H p = −Eq. The equations of motion of H = H 0 + H p are: with flow The dissipation function with initial distribution f 0 = e −βH 0 /Z 0 is given by: since the vector field is Hamiltonian and Λ vanishes. Moreover, we see that Ω and since p 0 = dq dp f 0 (q, p)p = 0, q 0 = dq dp f 0 (q, p)q = 0 and p 2 0 = dq dp f 0 (q, p)p 2 = β −1 . The response formula (78) gives Analogously, the unperturbed dynamics, given by E = 0 in (76), yields: and from (79) that must be multiplied by p before taking the ensemble averages in (78) and (79). This gives a linear combination of monomials p, pq, pq 2 , p 2 q, p 3 , whose average with respect to · 0 vanishes, and a non vanishing term p 2 E ω 2 − cos 2 ωt + sin 2 ωt + cos ωt Thus, recalling that p 2 0 = β −1 , we have: On the other hand, from (84) with E = 0, we obtain that p(O • S t 0 ) is a linear combination of q 2 p, qp 2 , p 3 . This implies that βEp O • S t−s 0 0 = 0 and, from (79), that qp ∼ t = 0, showing that the perturbed and unperturbed responses substantially differ in the case of O = q p, unless E is particularly small.

Formal Thermodynamics
Following Qian and collaborators [29], we now obtain various new relations for the dissipation function. To this end, suppose that the system admits an energy function H : M → R, H ∈ C 1 (R N , R). Let us introduce the set Γ h in phase space, defined as follows: and its boundary: The set Γ h includes all points in phase space whose energy is equal to or lower than the threshold h, and ∂Γ h is its boundary, constituted by the phase points whose energy equals the threshold value. Introduce the measure of such sets w.r.t. the time evolving probability measure: where f t (x) is a sufficiently well behaved probability density. Then, the function where χ h is the characteristic function of Γ h , plays the role of a probability density in Γ h . Using Reynolds Transport Theorem plus the regularity of the surface ∂Γ h , one finds: where dσ is the surface element area. This quantity can be interpreted as the measure w.r.t. the density f t (x) of the constant energy surface ∂Γ h . The time derivative along a trajectory of the energy function H(x) is given by: hence we may define a formal heat flux through the boundary ∂Γ h as: If in addition, we introduce the formal temperature as the inverse of the logarithmic derivative of the constant energy volume w.r.t. the energy parameter h: we can relate the ensuing formal entropy production rate to the average of the dissipation function concerning the set Γ h : This result is readily obtained as follows: where the equality is guaranteed by Divergence Theorem, ∇H/ ∇H being the unit vector everywhere normal to the surface ∂Γ h , and by the identity ∇ · Then Equations (3) and (5), yield: that leads to: Therefore, we can now write: The first addend within the brackets can be reworked using the response formula (8), where Ω f 0 restricted to Γ h is considered like any other phase function: (92) while for the second contribution, we have: Rearranging the two contributions, the first integral in Equation (91) becomes: and finally: Letting h grow in Equation (88), these conditional averages converge to the average on M, yielding: lim where (13) has been used. This is consistent with the fact that in this model no "energy" can flow out of M, or in M from outside. Note that such a reasoning actually does not refer to a physical flow: in the first place, it occurs in phase space, rather than in real space; furthermore the dynamics does not need to describe any physical phenomenon. Therefore, the thermodynamic interpretation of the above relations should be taken with a grain of salt [43,44], while it may be suggestive in the analysis of dynamical systems, analogously to the thermodynamic formalism [33].

The Case with Constant or Slowly-Varying Volume Contraction
Assume that the probability density f 0 is continuous with compact support M, and take V(x) = V 0 (x) + V 1 (x), with V 0 , V 1 sufficiently regular, Λ 0 (x) ≡ divV 0 (x) = −λ 0 < 0. Assume also that the system is Ω T-mixing. To compute Ω f 0 V t up to O( ), observe that: where The second integral can be computed explicitly taking f 0 (x) = 0 on ∂M: As the divergence of V 0 equals the constant −λ 0 , and M f 0 (x)dx = 1, we have: The integral in the previous equation is, by the Divergence Theorem, equal to and Averaging Λ one has instead: since to first order, our vector field has a constant phase space volumes contraction rate.

Example 1.
The simplest case is obtained when the initial distribution is uniform, i.e., f 0 (x) =constant, as if the "equilibrium" vector field vanishes. In that case the integrals in (101) vanish, and

Example 2.
Suppose the coordinates in the initial ensemble are independent, i.e., with g i being one dimensional density functions. Then, and Equation (101) can be rewritten as follows: In the special case of Gaussian marginals we have: The simplest example can be obtained with = 0 and V 0i (x) = −λ i x i , λ i > 0, and assuming µ k = 0. Then, observing that the i-th component of S τ 0 x, can be written as (S τ 0 x) i = e −λ i t x i , Equation (105) becomes: where λ 0 = ∑ n i=1 λ i . In addition, recalling that The previous result could have been obtained directly from Equation (96), noting that: Equation (107) shows that the average dissipation function Ω f 0 V t diverges as t → ∞, because the dynamics has the global attractor x = 0. Thus, Ω f 0 V tends to the non vanishing phase space contraction rate −λ 0 . That is not the case for > 0.

Simple Attractors and Small Dissipation
, modulo a set of zero Lebesgue measure. Let f 0 be the initial probability density and f t its evolution at time t. Recalling that the basins of attraction are pairwise disjoint, for any continuous phase function C we can write: Taking the limit, by bounded convergence we have: This shows that the asymptotic mean of C is the weighted mean of the values of C on the attracting fixed points, with weights given by the µ 0 -probability measure of the corresponding basins of attraction.
In other words, the probability measure dµ t (x) = f t (x) dx converges weakly in time to ν = ∑ k i=1 w i δ x i , with δ x i the atomic measure on x i . In this case it is easy to compute the correlation function of A and C w.r.t. µ 0 : that is used to define T-mixing and, taking C = Ω f 0 , to define ΩT-mixing, cf. Equations (10) and (11). The t → ∞ limit of (110) then yields: (111) and, assuming that also Ω f 0 is continuous, since Ω f 0 0 = 0.

Remark 1.
For k > 1, the response from a single initial condition is not unique and depends on the initial distribution f 0 . Hence, T-mixing does not hold, since it requires the long time averages to be independent of the initial condition [26]. On the other hand, ΩT-mixing holds, since it is necessary and sufficient for the convergence of the full phase space averages to their asymptotic limit, i.e., (109).
Now, consider a case with small dissipation: where is small, and divV 0 (x) = 0 ; Then, the unperturbed dynamicsẋ = V 0 (x) preserves the volume and the initial density, since Ω f 0 V 0 (x) ≡ 0, while the dissipation function only depends on the perturbation V 1 : Let the initial probability density f 0 be positive for every x in the phase space M, a condition at times called ergodic consistency [30]. In this case, the probability density at time t is given by: where S τ is the flow of (113). Taking α ∈ (0, 1) and t = 1/ α , it can be shown that the difference between the evolved density and the initial one is O( 1−α ), in the -dependent time interval bound by −t , provided Ω f 0 V is continuous in M. In particular, one has: which defines the maximum magnitude of Ω f 0 V 1 in M. Consequently, for 0 < t ≤ t : We can use this inequality to bound from above the variation of an observable of the system at hand. For a continuous function B, one obtains:

Conditions on Dynamics and Probability Density
As observed in [27], the dissipation function can be obtained starting from the Liouville equation: and factorizing f t (x) in the right hand side. The result is Equation (2) with Ω f t (x) defined as in (3). The standard conditions of physical interest, under which this operation can be performed are the following: • the vector field should be everywhere differentiable in the phase space M for Λ to exist; • the initial density f 0 should be everywhere positive in M for its logarithm to exist; • the initial density f 0 should be everywhere differentiable in M for its gradient to exist.
This also guarantees that the dissipation function be continuous, which is what is used to derive Equation (8) from Equation (7), via Equation (5). In the physical literature, these requirements are satisfied if (a) the condition sometimes called ergodic consistency holds [30]f 0 should be the equilibrium ensemble for the dynamics subjected to the same constraints of the nonequilbrium dynamics-and if (b) the particles interactions are sufficiently regular. The effects of discontinuities can be controlled if the dynamics are piecewise regular (as in billiards), or the probability decays sufficiently fast when approaching singularities. A full mathematical treatment of these aspects is far from trivial for physically relevant particle systems, and the success of the theory can only be judged a posteriori, on a case by case basis.
At the same time, various mathematical aspects can be understood considering simple examples, that can be explicitly solved. Let us focus on simple cases in which the above conditions are violated. This is important both to delineate the range of applicability of the theory based on the dissipation function and, at the same time, to find corrections where needed.
When f t (x) is not differentiable, one may replace (119) by its weak form, i.e., C ∞ 0 (M) being the set of smooth functions with compact support in M. The non-differentiability, of f t , however, is not the only issue. In what follows we take f 0 and f t with compact support K 0 and K t in M, while M is not necessarily compact.

Example 3. Consider the dynamical systemẋ
and the following initial density: that is compactly supported and non differentiable at 1 and −1. As the initial condition x evolves like S t x = e −t x, f t is given by Then, Equation (119) does not apply, while Equation (120) does. Indeed, the left hand side and the first term in the right hand side of (120) are The phase integral in the right hand side of (120) yields: and the corresponding time integral becomes: Then, integrating by parts and using Leibniz's formula in the last integral, we get: and Form (124) we get A = B + C, that proves (120).
The probability with density f t (x) converges weakly to the delta function at x = 0, since for any continuous O, the mean value theorem yields: for one x t ∈ (−e −t , e −t ). Can this result be obtained via the response formula (8)? Because Λ(x) = −1, and d f 0 /dx = 0, the dissipation function takes the form: showing that Ω In Ref. [27], the equality Ω f t f t = 0 was obtained from the conservation of probability and from the interchange of the time derivative with the integration in phase space: This exchange is legitimate under the joint continuity w.r.t.
x and t of f t and of ∂ f t /∂t, which is not the case of our example. Consequently, Equation (8) leads to the nonsensical expression: which diverges linearly with t, rather than converging to O(0), because O(e −t x) → O(0). This is related to the fact that Ω f 0 0 = 0, due to the non differentiability of f 0 . If, on the other hand, one has 1], one obtains: Then, taking for instance O(x) = x n , which can be used for analytic functions, one has: for even n, which correctly converges to O(0) = 0, and 0 for odd n, which equals O(0) at all times. The role of Equation (13) is further clarified noticing that for a finite t one can write: Then, assuming the t → ∞ limit and the phase integration can be exchanged in Equation (135), the previous equation yields: where O is the time average of O. This shows that Equation (9) requires Ω f 0 0 = 0, i.e., that either Ω f 0 vanishes or changes sign in M.
On the other hand, one realizes that this is not sufficient in general. For instance, despite the correct result (134), the blind application of Equation (5) for the case of Equation (121) with density 3(1 − x 2 )/4 in [−1, 1], yields an incorrect evolved density f t . In this case, we have: and which implies: that is not normalized on [−1, 1]. This difficulty is not related to the violation of Equation (13), but to the fact that the support of f t changes with time or, equivalently, that Equation (5) attributes a positive probability to points outside [−1, 1]. An analogous, even more macroscopic, situation arises if M = [1,2] and f 0 (x) = 2(3 − x)/3, which yields V(1) f 0 (1) = V(2) f 0 (2), hence yields (13) with non vanishing V(1), V(2) and f 0 (1), f 0 (2). This stresses that dynamical equations and phase space must be properly chosen. The point is that the preimages S −s x of part of the points in M, for some s > 0, lie outside M. In these cases, Equation (5) requires that only the contributions of the points x with S −s x ∈ M be included in the calculations, as discussed in Ref. [32].

Example 4. Consider the dynamical systemẋ
that has two asymptotically stable fixed points x − = −1 and x + = 1, and an unstable fixed point, x 0 = 0. Take the invariant set M = [−1, 1] as the phase space of the system. The flow is given by: Assuming f 0 (x) = 1/2, uniform on M, the dissipation function is expressed by: Ω f 0 (x) = −Λ(x) = 3x 2 − 1, hence it changes sign in M. The evolved distribution is given by: If, on the other hand, f 0 (x) = 3 4 (1 − x 2 ), we have that also changes sign, and the density at time t becomes: From (141) and (143) one can see that in both cases dµ t = f t (x)dx converges weakly to 1 2 (δ −1 + δ 1 ). This is correctly so because the flow vanishes at the boundaries of the phase space: M is invariant; hence the points inside M have no preimage outside M. In this case, the non-differentiability of the uniform density in R has no bearing on the theory.

Example 5.
Let us consider the one-dimensional dynamical system: and denote V 0 (x) = −x and V p (x) = x 3 . This system has three fixed points: x 0 = 0, which is asymptotically stable, and the unstable points x − = −1/ √ and x + = Then, expanding with respect to , This allows us to check that the probability measure dµ t (x) = f t (x; ) dx converges weakly to the delta function at x = 0 as t → ∞ The linear approximation for small can be obtained, in analogy with the Hamiltonian case, cf. (54), by truncating the exponential in (47): which gives:f It can also be observed that the linear approximation at = 0 of the exact evolved density (147), does not coincide with Equation (149), but it can be obtained from Equation (149) approximating ln x to first order in x.
Given that Ω , the exact response of ensemble averages is given by: that, on the basis of the previous remarks, converges to O(0), as t → ∞. As observed in Section 3, the linear approximation analogous to the GK response is obtained by replacing the exact flow S t x with the linearized one, S t 0 x = x e −t , in Equation (151). Now, let A and B be smooth functions, and f 0 a probability density with compact support K in M , containing x = 0 in its interior, for instance K = [−1, 1]. In order to check T-mixing, see (11), let us write: which follows from the uniform convergence of S t x to 0, from the boundedness of A and B in any compact set and from the continuity of A. This proves the T-mixing property (11) for the present example. The same reasoning shows that the system is also ΩT-mixing.

Example 6.
Consider the one-dimensional dynamical system: that is the time reversed version of (144). The point x 0 = 0 is now unstable, while x − = −1/ √ and where S t 0 x = xe t . The expansion at = 0 has an error that grows without limit as t increases, which shows the unsuitability of the linear response. In this case, always assuming f 0 being Thus the measure dµ t (x) = f t (x; ) dx converges weakly to 1 2 δ −1/ √ + 1 2 δ 1/ √ .

Example 7. Let us considerẋ
and assume the invariant interval M = [−1, 1] as the phase space. For = 0 the fixed points are x 0 = 0, x − = −1 and x + = 1. For > 0, the point 0 is unstable, while −1 and 1 are stable. The opposite holds for < 0. For = 0, all the points are fixed. The flow is given by: S t x = xe t 1 + x 2 (e 2 t − 1) .

Concluding Remarks
In this paper we have reviewed the developments that led from nonequilbrium molecular dynamics (NEMD) to the fluctuation relations (FR) and to the exact response theory, which develops the previous formalism based on the so-called Transient Time Correlation Function [6]. Within this framework, some authors have adopted the phase space contraction rate −Λ as a definition of the steady state entropy production rate σ [23,41], after −Λ and σ had been observed to be proportional to each other, in the original paper on fluctuations relations [36]. That paper dealt with the Gaussian thermostatted iso-energetic model of a shearing fluid known as SLLOD. In that case, the steady phase space contraction and entropy production rates could indeed be directly related to each other, something that in NEMD models had already been pointed earlier [19,51,52]. Mostly the interest had been on steady state state averages of such quantities. It was then realized that such an identification is justified only for steady state averages of certain NEMD models and that, in general, fluctuations require −Λ to be distinguished from the energy dissipation rate [24,25,53]. At the same time, the function that represents energy dissipation in molecular dynamics has been identified as the dissipation function Ω, the function that obeys the FR. For the transient FR this happens quite generally; for the steady state FR it is required that correlations w.r.t. the initial distribution do not grow too fast, something implied by the T-mixing condition [26].
Through a number of examples, including NEMD models as well as simple exactly solvable models, we have illustrated these facts and how the exact response formula, based on Ω, applies. In particular we have analyzed the relation between Ω and −Λ, noting that in the analysis of generic dynamical systems the two quantities are by and large equivalent, while in the case of particle systems the first corresponds to a real observable. We have then compared the exact response formula with the one of the Green-Kubo linear theory, quantifying their difference, that typically grows with the size of the perturbation and of time, but also showing that in some simple case they yield the same result.
We have then investigated "formal thermodynamics" along the lines of Ref. [29], introducing a formal temperature and a formal energy flow, and using Ω instead of −Λ as a notion of formal entropy production rate. Analogously to earlier works, e.g., Refs. [54][55][56], this means that the phase space is treated as real space, and its different regions as different locations in real space. Clearly, the 2dN-dimensional phase space of an N particle system in d dimensions and the phase points that flow in it are totally different from the 2d-dimensional real space and the particles that flow through it [43,44]. (A point in phase space has no extension nor "inertia", it does not interact with other points, and it equally represents the microstate of a large or small system. Any other point represents another system that can be arbitrarily close or far from it. A particle in real space has a finite extension, interacts with nearby particles, and follows Newton's laws with its inertia. Just this has numerous consequences. For instance, while the Gibbs entropy represents thermodynamic entropy only at equilibrium, the Boltzmann entropy works also in time dependent situations, Refs. [44,57]). Therefore, in general, a direct identification of flows and related quantities in phase space with their physical counterparts is not possible. Nevertheless, this analysis is reminiscent of the thermodynamic formalism, which is important for dynamical systems in general [33]. Furthermore, it may find direct physical application in systems of non-interacting particles, for which the factorization of the phase space distributions makes phase space and real space almost the same. Such systems may not behave thermodynamically, but they are of greater and greater interest in e.g., bio-and nano-technology (cf. transport in highly confining media), for photons in optical wave-guides or random media, for cold atoms etc.
Finally, we have observed that the exact response formula (8), mainly rests on the regularity of the vector fields, and on the ergodic consistency of the initial distribution, which ensure the regularity of the dissipation function as well. These mathematical conditions need further scrutiny, in more realistic cases, but are relatively minimal, hence hint at a wide applicability of the theory. The investigation of physical applications of this theory is just at the beginning; many are indeed the cases in which perturbations cannot be taken as small compared to the reference signals; they range from nano-technology to geophysical and even astrophysical phenomena [58].
Author Contributions: All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication. All authors have read and agreed to the published version of the manuscript.