Using the fluctuation-dissipation theorem for nonconservative forces

An equilibrium system which is perturbed by an external potential relaxes to a new equilibrium state, a process obeying the fluctuation-dissipation theorem. In contrast, perturbing by nonconservative forces yields a nonequilibrium steady state, and the fluctuation-dissipation theorem can in general not be applied. Here we exploit a freedom inherent to linear response theory: Force fields which perform work that does not couple statistically to the considered observable can be added without changing the response. Using this freedom, we demonstrate that the fluctuation-dissipation theorem can be applied for certain nonconservative forces. We discuss the case of a nonconservative force field linear in particle coordinates, where the mentioned freedom can be formulated in terms of symmetries. In particular, for the case of shear, this yields a response formula, which we find advantageous over the known Green-Kubo relation in terms of statistical accuracy.

The linear response of a classical equilibrium system to a potential perturbation U ptb applied for time t > 0 is given by the fluctuation-dissipation theorem (FDT) [1][2][3][4], where A is an observable of interest, k B is Boltzmann's constant, T is temperature, and · · · ptb and · · · indicate averages in the perturbed and equilibrium system, respectively. The stationary limit of formula (1) can be derived from the equilibrium distribution the system relaxes to.
In this work, we exploit a simple method for computing the linear response to a nonpotential perturbation via FDT for certain observables, using the freedom of adding forces whose work does not couple to the considered observable. Note that a similar freedom has been discussed in Ref. [22]. We illustrate that, for a force perturbation linear in coordinates, the mentioned freedom can be formulated in terms of symmetries. We study in detail the case of a Brownian system perturbed by simple shear flow, finding a response formula [Eq. (11) below], which is an alternative to the classical Green-Kubo relation [Eq. (12) below]. Using numerical simulations, the formula is found to have a lower variance for all cases studied, making it advantageous in terms of statistics. Consider a classical system of N interacting particles, subject to external potentials and coupled to a heat bath at temperature T , in thermal equilibrium at time t = 0. For time t > 0, the system is perturbed by nonconservative forces {F ptb i }, with F ptb i acting on particle i at position r i . Because the equilibrium state is time symmetric, the linear response of A is related to the work done on the system [3,23], (2) If F ptb i is a conservative force, i.e., F ptb depends only on the states, and Eq. (1) follows from Eq. (2).
Notably, response relation (2) displays a freedom in F ptb i when computing the perturbed A: It allows adding perturbation forces G ptb i whose work does not couple to the observable A, i.e., without changing the response of A. Thus, if a force G ptb i obeying Eq. (3) exists such that adding the two forces results in a potential U ptb , then, according to Eq. (2), the response of A to the nonconservative force F ptb i is equivalent to the response to the potential U ptb and given by formula (1). Exploring this possibility of restoring an FDT is the content of this paper. We investigate the specific case of a force field linear in r i , with the tensor κ independent of particle positions. If κ is symmetric, F ptb i derives from a generalized harmonic potential. The case of interest is that κ is not symmetric, such that F ptb i of Eq. (5) is not conservative. One natural way of exploring the above-mentioned freedom is by using the transpose of κ, i.e., it is promising to use The sum of F ptb i and G ptb i is then immediately found, where the potential is identified as How to satisfy Eq. (3)? Many cases that do so can be identified on the basis of symmetries, as we demonstrate by regarding κ = κ 12x ⊗ŷ (withx,ŷ, and ⊗ denoting unit vectors and the tensor product, respectively), i.e., shear forces (see Fig. 1 for an illustration) [24]. From Eq. (8), the corresponding potential reads being a potential with one stable and one unstable direction in the xy plane (see Fig. 1). One direct way of fulfilling Eq. (3) is restricting to systems and observables which are symmetric under interchange of the x and y coordinates. These are systems for which interaction and external potentials remain the same under interchange {x i } ↔ {y i }, and observables which remain the same under interchange {x i } ↔ {y i } and {v ix } ↔ {v iy } (where v ix denotes the x component of the velocity of particle i). Then condition (3) is fulfilled by symmetry [25]. For example, spherically symmetric potentials and observables or the xy component of the stress tensor [4] comprise these symmetries. Substituting Eq. (9) into Eq. (1), we find that, for these cases, the linear response to shear forcing is given by Formula (10) thus provides the response to a nonconservative force (shear force) via FDT. Many models treat forces and external flow driving in an identical manner. Regarding Brownian particles with mobility µ [see Eq. (13) below], an external flow velocity field V(r) gives rise to a force F ptb , neglecting hydrodynamic interactions. The mentioned shear forces then translate to shear flow by identifying κ 12 =γ µ (with shear rateγ), and, under the above symmetries, we obtain for the shear perturbed A(t) (γ) , Response relations (11) and (10) are our main results. We note that Eq. (11) has been derived in Ref. [17] for a single overdamped Brownian particle. A prestage version of Eq. (11) has been given in Ref. [27], and we discuss the relation at the end of this paper. Equation (11) is an alternative to the classical Green-Kubo relation for shear which, for the case of overdamped Brownian particles, reads as [see Refs. [2-4, 7, 18, 28, 29] for various Green-Kubo relations, and Ref. [18] for formula (12) in particular] where σ xy is the xy component of the stress tensor defined as ix being interaction and external forces, respectively, acting on particle i in direction x. One advantage of Eq. (11) over Eq. (12) is the absence of a time integral. Another is that forces do not have to be measured.
Seeking a numerical example, we turn to interacting overdamped Brownian particles in two space dimensions, following Langevin dynamics [3,23], where µκ · r i is the shear velocity, imposed at t > 0, with the shear-rate tensor µκ =γx ⊗ŷ.
are interaction forces, chosen to arise from a screened Coulomb potential, with interparticle distance r ij ≡ |r i − r j |, coupling strength Γ, and interaction range r c . The external force follows from a harmonic potential, with spring constant k. f i is a Gaussian white noise, where I is the identity matrix. We set k B T = r c = µ = 1, and Γ = 25. N , k,γ, and the number of independent noise realizations C for performing averages are varied between measurements. The dynamics is simulated using the Euler method. We choose A = N i=1 x i y i , which is the lowest nontrivial moment of the particle distribution.
Since the system and A are xy symmetric, condition (3) is fulfilled and formula (11) is valid.
We compute A(t) (γ) − A via three different routes: by (i) applying finite shear, (ii) using equilibrium correlations according to the Green-Kubo formula (12), and (iii) using equilibrium correlations according to Eq. (11) (labeled "FDT" in the figures). Figure 2 compares these as a function of time t after start of shear. For small shear rate (main plot), all methods agree, thereby confirming formula (11). For large shear rate (inset plot), the deviation from the linear response is evident, also regarding the form of the response curve, which shows a characteristic "overshoot," i.e., a nonmonotonic behavior as a function of time, which has also been observed in sheared bulk systems [30]. Snapshots for equilibrium (black particles) and sheared (orange particles) systems illustrate the change of shape of the cluster from circular to ellipsoidal: A = 0, but A(t) (γ) ≥ 0. Panels (a) and (c) of Fig. 3 show the dependence on the confinement strength k and the number of particles N of the steady-state response, again confirming agreement between the three methods. From fits to the data, the response follows the scaling ∝ k −1.48 (compared to ∝ k −2 , obtained analytically for N = 1) and ∝ N 1.55 (for N 4).  Fig. 3 show the corresponding variance, related to the statistical error of a single measurement using the different methods [31]. It shows a notable difference between the methods following scaling behaviors of ∝ k −0.64 N 0.82 , ∝ k −1.48 N 1.54 , and ∝ k −1.50 N 1.54 , respectively. The Green-Kubo relation and formula (11) scale similarly, but the latter has a notably lower variance. Table I compares scaling behaviors of the relative variance (variance divided by the mean) for the three methods. The relative variance of the directly measured response grows with k and decreases with N . For the Green-Kubo relation (12) and for Eq. (11), the relative variance hardly depends on k and N , indicating that the statistical efficiency of Eqs. (12) and (11) is invariant with respect to changes of the effective system size and density, highlighting an interesting property of the linear response approach. For the set of parameters we used in our simulations, Eq. (11) has the lowest variance. Comparing it to the Green-Kubo relation (12), it thus needs a much smaller number of independent runs (roughly a factor of 100 here, estimated from the variance and the  central limit theorem), which, additionally to the mentioned absence of integration, is advantageous.

Panels (b) and (d) of
Finally, the Langevin equation (13) allows one to give more insights into the nature of Eqs. (11) and (12). Expanding the corresponding path action in shear rateγ yields for the linear response [27] The second term on the right-hand side of Eq. (17), containing the stress tensor, stems from the time-symmetric part of the expanded action, and yields (the half of) Eq. (12). The term containingẋ i y i is time antisymmetric, and yields, after adding the transpose shear field, Eq. (11). Because the equilibrium state is time symmetric, the two terms in Eq. (17) are identical [32]. This discussion finally highlights another advantage of FDT and Eq. (10): The form of Eq. (10), being based on the time-antisymmetric part, related to the above-mentioned work, is system independent, while Eq. (12) takes different forms in different systems [33].
Using the freedom of adding forces whose work does not couple to the considered observable, we found that the linear response to nonconservative forces can be computed from FDT. Compared to standard approaches, application of this concept typically involves simpler quantities to be measured [e.g., positions in Eq. (11) versus forces in Eq. (12)], does not require time integrals, and necessitates a smaller number of independent measurements. This is expected to make it advantageous for both simulations as well as experiments. Future work may address bulk systems via Eq. (11) in order to obtain the shear viscosity and to connect to "Einstein relations" for such viscosity [4,[34][35][36][37][38][39].
We greatly acknowledge collaboration with Christian M. Rohwer and Alexandre P. Solon during the initial stage of this project. We also thank Rosalind J. Allen, Thomas Voigtmann, and Patrick B. Warren for valuable discussions, and thank S. Dietrich and Alexandre P. Solon for comments on the manuscript. K.A. is supported by Studienstiftung des deutschen Volkes and