Relaxation of a system to equilibrium is as ubiquitous, essential, and as poorly quantified as any phenomena in physics. For over a century, the most precise description of relaxation has been Boltzmann's H-theorem, predicting that a uniform ideal gas will relax monotonically. Recently, the relaxation theorem has shown that the approach to equilibrium can be quantified in terms of the dissipation function first defined in the proof of the Evans-Searles fluctuation theorem. Here, we provide the first demonstration of the relaxation theorem through simulation of a simple fluid system that generates a non-monotonic relaxation to equilibrium.

The behaviour of objects in our everyday world, be it coffee cooling or gases expanding, suggests that a system will naturally relax to equilibrium given sufficient time. This is critically important, as once a system is in equilibrium, we can define important thermodynamic properties such as temperature, entropy, and free energy. Indeed this seems a natural consequence of the second law of thermodynamics. In spite of this, quantifying how a system relaxes to equilibrium is far from simple.

Boltzmann's H-theorem was the first, and for a long time the only quantitative description of relaxation. Starting with a spatially uniform ideal gas and assuming molecular chaos, it can be shown that there is a function H which will decrease with time,^{1}

Physically, the restriction on the sign of the H function means that any relaxation must be monotonic. The H function is connected to the Gibbs entropy (S) of an ideal gas at equilibrium, (S = −*k*_{B}H), and it can be shown that *d*H/*dt* = 0 when the system has a Maxwell-Boltzmann distribution suggesting that the equilibrium distribution for the ideal gas is the distribution of maximum entropy for the system, and that this distribution is stable with respect to time. This relation emphasizes the connection between the relaxation and the second law of thermodynamics. Although this approach can be extended to other systems that feature molecular chaos, applying the H-theorem without this property has been extremely difficult.^{2,3}

Modern developments in understanding relaxation have largely been based on mathematical physics. It is possible to show that if a system is mixing then it will relax to a unique ergodic state.^{4} In practice, however, very few systems can be proven to have these properties, though they can generally be assumed.^{5} More of a problem is that without being able to prove these properties, this approach has little to say about *how* the system relaxes to equilibrium.

Ultimately, the difficulty is that the approach to equilibrium is inherently a non-equilibrium process, and there are very few analytic expressions in thermodynamics that apply outside of equilibrium or local equilibrium regimes. In the last two decades, there have been major advances in non-equilibrium thermodynamics with the development of the fluctuation relations: exact expressions that apply to finite size, non-equilibrium systems.^{6} At the core of these advances has been the development of a new response function, the dissipation function (Ω_{t}), that has a pivotal role in describing and quantifying the non-equilibrium behaviour of thermodynamic systems.

The dissipation function is the argument of a number of new relations such as the fluctuation theorem, the non-equilibrium partition identity, and the dissipation theorem.^{7–9} The time integral of the dissipation function (total dissipation) is defined as the logarithm of the probability ratio of observing two sets of trajectories originating in the initial distribution that are time reverses of each other:

and the instantaneous dissipation function is therefore

Here, **Γ** ≡ {**q**_{1}, **p**_{1}, …**q**_{N}, **p**_{N}} is the phase space vector of the system, Ω_{t}(**Γ**(0)) is the total dissipation for a trajectory originating at **Γ**(0) and evolving for a time *t*, Ω(**Γ**(*t*)) is the instantaneous dissipation function, *f*(**Γ**(0), 0)*d***Γ**(0) is the probability of observing a system in an infinitesimal region around **Γ**(0) in the initial system distribution, and **Γ***(*t*) is the result of applying a time reversal map to **Γ**(*t*). The dissipation function (Ω_{t}) generally behaves like an entropy function: it is positive when the system is exhibiting behaviour consistent with the second law of thermodynamics, and it is negative when the system is exhibiting unusual anti-second law behaviour. It is however a path function, and not a state function. This means that to study overall system behaviour we need to study either the distribution of the dissipation function such as in the Evans-Searles fluctuation theorem (ESFT),^{7}

or the ensemble average of the dissipation function. Note *P*(Ω_{t} = *A* ± *dA*) is the probability of observing a trajectory with a dissipation value infinitesimally close to *A*. Importantly, the ensemble average of the total dissipation function obeys a second-law-like expression, 〈Ω_{t}(**Γ**(0))〉 ⩾ 0 , ∀*t* > 0, so that while any individual system trajectory may be anti-dissipative, the ensemble behaviour is not.^{10}

The relaxation theorem quantifies the relaxation of a system to equilibrium through the dissipation function, rather than the H-function. To derive the relaxation theorem, the system must have time decay of correlations, an initial distribution that is even with respect to time reversal, and exhibit ergodic consistency, that is, the dynamics must constrain the occupied volume of phase space with time to a subset of the initial phase space volume.^{11,12} From here, several applications of both the second law inequality and the dissipation theorem can be used to derive the relaxation theorem. The dissipation theorem is a restatement of the transient time correlation function for the system with the dissipation function as a core argument,^{9}

where *B*(**Γ**) is a phase function, Ω(**Γ**(0)) is the instantaneous dissipation at time 0, and 〈…〉 is an ensemble average, that is, an average over all the available points in the initial non-equilibrium distribution.

Like the H-theorem, the relaxation theorem describes the average behaviour of systems that satisfy a set of conditions. For a system with these properties, the relaxation theorem predicts that:^{11}

The instantaneous dissipation will be zero for a system in equilibrium (Ω(

**Γ**(*t*)) = 0).The state where the instantaneous dissipation remains 0 with time is unique and is the canonical distribution for a thermostatted system.

The average instantaneous dissipation will go to zero as the system approaches equilibrium (lim

_{t → ∞}〈Ω(**Γ**(*t*))〉 = 0).If the system relaxes conformally,

^{13}the approach to equilibrium will be monotonic (〈Ω(**Γ**(*t*))〉 ⩾ 0 ∀*t*> 0).For non-conformal relaxation, the

*average*total dissipation will still always be greater or equal to zero (second law inequality: 〈Ω_{t}(**Γ**(0))〉 ⩾ 0 ∀*t*> 0).

In other words, as a system relaxes the average instantaneous dissipation will go to zero and the total dissipation will be larger than zero. This does not prescribe monotonic relaxation such as the H-theorem: there is no constraint that the total dissipation must always increase. To follow the relaxation of a system, we only need to measure the instantaneous dissipation with time.

A very simple system for testing the various fluctuation relations is based on an optical trapping experiment called the capture experiment.^{14} In this system, a single particle is bound to a harmonic potential created by a laser, which is equivalent to tethering the particle to a point in space with a Hookean spring. The particles in the fluid constantly interact with the bound particle causing it to move in response, while the harmonic force prevents the particle from moving too far from the trap centre. These two competing forces mean that at equilibrium, the trapped particle will have a Gaussian distribution of positions and velocities, defined by the strength of the trap and the temperature of the fluid. In our system, we model this with a two-dimensional fluid of Weeks-Chandler-Anderson particles surrounded by a double layer box of harmonically bound, Nose-Hoover thermostatted particles further bounded by periodic boundary conditions.^{15,16} In the centre of the box, one of the fluid particles is bound by the harmonic trap that drives the experiment.

To perturb the system from equilibrium we discontinuously change the strength of the trap at *t* = 0^{+} from *k*_{0} to *k*_{t}, and allow the system to relax. This means the system starts in a known distribution with *k*_{0} and relaxes to the known distribution of *k*_{t}. The equations of motion for the system are

where *i* is the particle index, *m* is the mass of the particles, δ is the Kronecker delta, *S*_{w} is the thermostat switch that is 1 for the wall particles and 0 for the fluid and trapped particles, **F**_{I, i} is the intermolecular force acting on the particle, **F**_{w, i} is the harmonic force constraining the wall particles to their positions, *Q* is the thermal mass of the thermostat, *T* is the thermostat target temperature, and *T*_{k} is the kinetic temperature of the walls. The Nose-Hoover thermostat is an integral feedback thermostat applied to the wall particles. While it is a synthetic thermostat, by placing it in the walls of the system it acts indirectly on the fluid system in the same manner as a real thermostat.^{17}

The system begins in equilibrium with the thermostat and a Hookean spring of trapping constant *k*_{0}. From the initial distribution function and the equations of motion, Eqs. (6)–(8), we can derive our dissipation function,^{14}

where β = 1/*k*_{B}*T*, **q**_{1} is the position of the trapped particle relative to the harmonic trap centre, and **p**_{1} is the moment of the trapped particle. From the second law inequality, we know that the average total dissipation will be positive. Therefore, given the relative values of *k*_{0} and *k*_{t} we can predict whether

*k*

_{0}= 2 and

*k*

_{t}= 8, and therefore expect that the particle will move closer to the trap on average.

^{18}

If we examine the instantaneous dissipation function in this case, we expect a positive value when the particle is moving towards the trap centre, and a negative value when it is moving away. In order for the system to be in equilibrium, the probability of these two states must be equal, and therefore the average of the dissipation function would be zero as predicted by the relaxation theorem. If the trapped particle interacted with an idealised thermal bath that uniformly dissipated the extra energy until the trapped particle was back in thermal equilibrium with its surroundings, we would expect to observe monotonic relaxation. In our system, the fluid provides energy storage modes in the kinetic and potential energy of the surrounding fluid and the walls, corresponding to local heating, that can be returned to the trapped particle instead of being dissipated in the thermostat. Similarly, the Nose-Hoover thermostat can overshoot equilibrium if shocked, and then put heat back into the system. This can lead to non-monotonic relaxation in the system, which we aim to quantify with the dissipation function.

Before we test the relaxation theorem, we first demonstrate that the underlying dissipation relations (the Evans-Searles fluctuation theorem and the dissipation theorem) are obeyed for this system. In Figure 1(a), we plot the ESFT for the system at the end of the simulation, and see extremely good agreement with theory. In Figures 1(b) and 1(c), we take advantage of the fact that dissipation is a phase function to test the dissipation theorem. In Figure 1(c), we plot the difference between the two functions and the combined standard error for the two functions: if the two functions are in agreement we would expect the difference to remain between the error lines for the majority of points, which it does. This is the expected result as the validity of the dissipation theorem is necessary for the relaxation theorem to hold.

Finally, we plot the ensemble average of the total dissipation function against time. The system relaxes towards zero instantaneous dissipation at long times. This agrees with the prediction of the relaxation theorem that as the system goes to equilibrium, the average of the dissipation function must also go to zero. What is interesting about this system is that the dissipation overshoots, that is, for periods of time the average instantaneous behaviour of the system is anti-dissipative. Despite these regions, the average of the total dissipation function is always positive definite, and the value converges at long time.

From these results we can see that the Relaxation theorem combined with the second law inequality accurately models the non-monotonic approach to equilibrium. Furthermore, it can be shown that the dissipation function defined from any time *t*_{1}, (Ω_{τ}(**Γ**(*t*_{1}), *t*_{1})), where the second *t*_{1} represents the time at which the phase space distributions in Eq. (2) are evaluated and τ is the length of the trajectory, can be connected to the dissipation function defined from *t* = 0, (

_{τ}(

**Γ**(

*t*

_{1}),

*t*

_{1}) ⩾ 0).

^{19}This means that if you defined the dissipation with respect to the first non-equilibrium peak in the dissipation function, approximately

*t*

_{1}= 0.256, the total dissipation function would still be positive.

With this successful demonstration of the relaxation theorem, it is appropriate to mention a less obvious consequence of the theorem. For a system that has an initial distribution that is even in the momenta, and whose dynamics are T-mixing (i.e., in the infinite time limit the ensemble average of any smooth phase function is constant), the final equilibrium distribution is ergodic. This means that for deterministic systems, the phase space distribution is unique, and there is no room for alternatives to the canonical distribution as has occasionally been speculated.

The relaxation theorem represents a significant advance in understanding the relaxation to equilibrium. It provides a theoretical description of the relaxation to the Canonical distribution. Practically, the dissipation function is generally both calculable and measurable in experiments and simulations. From the form of the dissipation function we can make predictions about the direction of relaxation and the convergence of the average dissipation towards zero can be used as a measure of the degree of relaxation of the system.

## REFERENCES

*f*(

**Γ**,

*t*) = exp ( − β

*H*+ λ(

*t*)

*g*(

**Γ**))/

*Z*, ∀

*t*) and the deviation function,

*g*, is a constant over the relaxation.

*t*= 1 × 10

^{−3}, T = 1, ρ = 0.3, and 100 000 trajectories.