Extrapolation to nonequilibrium from coarse grained response theory

Understanding properties of nonequilibrium systems is an ambitious goal of modern statistical physics Seifert (2012), and here, the fluctuation dissipation theorem (FDT) is of fundamental importance: It relates the linear response of a system to its thermal fluctuations in the equilibrium state Kubo and Tomita (1954); Kubo et al. (2012). This insight is of practical benefit in solid state physics Phillips (2012) as well as in classical systems.

The FDT holds close to equilibrium, and extending it to far from equilibrium has been the subject of intense research. The case of small perturbations of far-from-equilibrium states has been analyzed in various works Seifert (2012); Ruelle (1998); Harada and Sasa (2005); Speck and Seifert (2006); Blickle et al. (2007); Chetrite et al. (2008); Marconi et al. (2008); Baiesi et al. (2009); Prost et al. (2009); Krüger and Fuchs (2009); Lippiello et al. (2014). Another direction aims at finding the nonlinear response, i.e., the response to strong perturbations. The derived formulas relate response functions to nonequilibrium correlation functions Yamada and Kawasaki (1967); Evans and Morriss (1988); Fuchs and M. E. Cates (2005), or to (higher order) correlation functions evaluated in equilibrium Kubo and Tomita (1954); Semerjian et al. (2004); Bouchaud and Biroli (2005); Andrieux and Gaspard (2007); Lippiello et al. (2008a); Colangeli et al. (2011); Lucarini and Colangeli (2012); Basu et al. (2015). The latter concept has been applied experimentally only recently Helden et al. (2016), where the second order response was obtained from an equilibrium measurement.

Extensions of FDT to far from equilibrium cases are typically plagued by a property, which is deeply inherent to nonequilibrium physics: Their application requires information about the interactions and dynamics of the system, so that in principle all degrees of freedom (or their nonequilibrium-distributions) have to be tracked during the measurement (see discussions in Refs. Basu et al. (2015); Seifert (2012)). This statement may be exemplified for colloidal particles, investigated in Ref. Helden et al. (2016): To apply second order response theory, the interaction potential of the particles and their dynamical laws have to be known (and monitored). It is this aspect of nonequilibrium response theory (the dynamical details mentioned in the abstract) which often restricts its applicability to systems with small number of degrees of freedom, and has prevented application to many body systems.

A general route for many body systems identifies a relevant subset of important (slow) degrees of freedom, and less relevant (fast) degrees are integrated out. Examples are the so called Mori-Zwanzig projection formalism Zwanzig (2001, 1960, 1961); Mori (1958, 1965) or Fokker-Planck- or Langevin equations Risken (1989); Dhont (1996). When applying such approaches to nonequilibrium cases, the integrated degrees of freedom are typically assumed to be in equilibrium.

In this Letter, we derive a response scheme which overcomes these issues: Starting from the microscopic description, we derive a nonlinear response relation for a small subset of coarse grained degrees of freedom, which is then used in an extrapolation scheme: Measurements near equilibrium, i.e., linear in perturbation, are used to predict responses further away from equilibrium, i.e., to second order in perturbation. The microscopic degrees do neither have to be tracked, nor are they assumed to be equilibrating fast, so that this scheme is applicable to many body systems. We demonstrate applicability in an exactly solvable jump process and in computer simulations of the 2D Ising model.

Coarse grained nonlinear response theory from path integrals – We consider a classical many body system which is in weak contact with an equilibrium thermal bath. Considering for example the Ising model (see below), nonlinear response theory, as e.g. given in Ref. Basu et al. (2015); Lippiello et al. (2008b), can only be applied if the Hamiltonian (e.g., nearest- or next to nearest neighbor interactions) and the dynamics (e.g., specific spin flip rules) are known, and if all degrees are tracked. Our goal is development of a nonlinear response method which can be applied by tracking a small number of degress of freedom, e.g., the order parameter in the Ising model, not necessitating knowledge about the details of the system.

To this end, we introduce a coarse grained description in terms of $n$ (experimentally trackable) macrostates, each containing several, uniquely assigned microstates. At any time $t$, the system is thus characterized by a unique macrovariable $X_{t}=0,1,\dots,n-1$ (e.g., the sign of the magnetization in the Ising model is described by two macrostates, $X_{t}=0,1$) which represents the coarse grained phase space. In the absence of perturbations, the system is in thermal equilibrium, and thus satisfies detailed balance and time reversal symmetry.

We aim to compute the response of the system to a perturbation, whose strength is quantified by the dimensionless parameter $\varepsilon$. The perturbation can for example be a force, an external field, or a change in the transition rates of a jump process. We restrict here to perturbations which are switched on at time $t=0$, but are otherwise time independent. We build on path integrals, in terms of which response theory has been worked out for the microsystem Colangeli et al. (2011); Basu et al. (2015); Helden et al. (2016): The probability weight $p(\omega)$ of a microscopic path $\omega$ in the perturbed process differs from its equilibrium weight $p_{\rm eq}(\omega)$. This is captured by the action $\mathcal{a}(\omega)$, i.e., $p(\omega)=e^{-\mathcal{a}(\omega)}p_{\text{eq}}(\omega)$. $\mathcal{a},$ which vanishes for $\varepsilon=0$, is expanded in powers of $\varepsilon$,

where the primes denote derivatives w.r.t. $\varepsilon.$ In the spirit of Refs. Baiesi et al. (2009); Colangeli et al. (2011); Basu et al. (2015), $\mathcal{a}=\mathcal{d}-\frac{1}{2}\mathcal{s}$ is split into a part symmetric under time reversal, $\mathcal{d}$, and an antisymmetric part $\mathcal{s}.$ We take the perturbation to be such that $\mathcal{s}$ is linear in $\varepsilon$, so that $\mathcal{s}^{\prime\prime}$ and higher derivatives vanish, which is a generic and useful case Kubo (1966); Baiesi et al. (2009); Seifert (2012). This may also be interpreted as a definition of the order of perturbation: $\varepsilon$ is the quantity, in which $\mathcal{s}$ is linear. For a perturbation via potential forces this means that the perturbation Hamiltonian is linear in $\varepsilon$.

The response of an observable, up to second order, can then be expressed in terms of equilibrium correlation functions involving combinations of $\mathcal{s}^{\prime}$ and $\mathcal{d}^{\prime}$ Baiesi et al. (2009); Colangeli et al. (2011); Basu et al. (2015) (we will refer to the corresponding response formula when introducing Eq. (9) below).

On the coarse grained level the probability $P_{ij}$ of the macro-path which connects $X=i$ at $t=0$ and $X=j$ at time $t$ ¹¹1The paths are sufficiently characterized by the parameters time $t$ and the initial and final states., follows from integration over microstates, and the corresponding macro-action $\mathcal{A}_{ij}$ is (in the following, we omit the time arguments for brevity, keeping in mind that, e.g., $\mathcal{A}_{ij}=\mathcal{A}_{ij}(t)$)

Here, $\int_{ij}$ denotes integration over only those micro paths $\omega$ which connect the macrostates $i$ (at $t=0$) and $j$ (at time $t$). Using the definition, $\int_{ij}d\omega p_{\text{eq}}(\omega)=P_{ij}^{\text{eq}}$, the right hand side of Eq. (2) may be expanded in a series of $\varepsilon$, to obtain the macroscopic analog of Eq. (1). For $\varepsilon=0$, the argument of the log is unity, and we use its expansion around that value to obtain, with the notation $\mathcal{A}=\mathcal{D}-\frac{1}{2}\mathcal{S},$

The first derivatives, $\mathcal{S}^{\prime}$ and $\mathcal{D}^{\prime}$, are thus given in terms of the microscopic counterparts, while notably, the coarse graining in general generates a finite $\mathcal{S}^{\prime\prime}$ in the last line of Eq. (3), although the microscopic counterpart $\mathcal{s}^{\prime\prime}$ is zero.

The expected value of a macro-observable $O(X)$ at time $t$ under the perturbation is given by the average over the macroscopic paths. Expanding $\mathcal{A}$ in powers of $\varepsilon$, we obtain, up to second order in $\varepsilon$,

Here $\langle\cdots\rangle$ and $\langle\cdots\rangle^{\rm eq}$ denote averages over the perturbed and equilibrium processes, respectively. Other terms in this expansion disappear because of time reversal symmetry of the equilibrium process, manifest here in the symmetry of the matrix $P_{ij}^{\rm eq}$. The last term in Eq. (LABEL:eq:soo) is not present in the microscopic version Basu et al. (2015), and it appears here because of the non-vanishing $\mathcal{S}^{\prime\prime}$ in Eq. (3). The extrapolation scheme proposed below is applicable if the last term in Eq. (LABEL:eq:soo) vanishes. In particular, it is instructive to consider perturbations which couple to the coarse grained variable $X$. One example is a perturbation potential $\varepsilon V(X)$, i.e., a perturbation potential which is sensitive to the macrostates. In that case, $\mathcal{s}^{\prime}(\omega)=\beta[V(X_{0})-V(X_{t})]$ Baiesi et al. (2009), with inverse thermal energy $\beta=(k_{B}T)^{-1}$. It is thus equal for all the micropaths connecting macro states $i$ and $j.$ Consequently, the term in the last line of Eq. (3) simplifies to

It immediately follows that $\mathcal{S}^{\prime\prime}=0$ in Eq. (3), and therefore, Eq. (LABEL:eq:soo) simplifies to a form

Equation (9), an intermediate result, is the response formula for the coarse grained phase space $X$. It is reminiscent of the microscopic version Basu et al. (2015), however here we obtained it for the croase grained variables. The left hand side is the nonequilibrium average of $O(X_{t})$, while the right hand side is an explicit expression in terms of the time dependent matrices $\mathcal{\mathcal{S}}^{\prime}$, $\mathcal{\mathcal{D}}^{\prime}$ and $P^{\text{eq}}.$ Important for this work is the interpretation of Eq. (9): It is worth appreciating that the second order response, given by the last term of Eq. (9), involves $\mathcal{S}^{\prime}$ and $\mathcal{D}^{\prime}$, which are the changes of these matrices to linear order in $\varepsilon$. This leads to the main result of the paper: Measuring the linear response of the system, i.e., measuring $\mathcal{\mathcal{S}}^{\prime}$ and $\mathcal{\mathcal{D}}^{\prime}$, is sufficient to predict the second order response from Eq. (9).

This extrapolation scheme does neither rely on the knowledge or tracking of integrated degrees of freedom, nor are they assumed to equilibrate fast (in contrast to Zwanzig-Mori approaches), and is thus applicable to many body systems with the caveat that the linear response needs to be measured. We illustrate the scheme in two examples.

Four state jump process – Let four micro-states, $A,\dots,D$ be connected with given jump rates, see sketch in Fig. 1. The coarse grained macrostates combine $A,B$ ($X=0$) and $C,D$ ($X=1$), respectively, so that $X$ is the phase space of a two state system ($n=2$), with $\langle X\rangle^{\rm eq}=\frac{1}{2}$ because of symmetry.

At time $t=0$, the system is perturbed by switching the forward rate of the central link from 1 to $e^{\varepsilon}$, while all other rates are left unchanged (see sketch in Fig. 1). Because we perturb the link connecting the macrostates, Eq. (9) can be used. We aim to find the responses up to the second order,

The response formula, Eq. (9), yields the predicted responses $\chi^{\text{rf}},$

Evaluating Eq. (11) in the extrapolation scheme, ${\cal S}^{\prime}$ and ${\cal D}^{\prime}$ need to be known. Therefore, the path weight $P_{ij}(t)$ is measured in linear response (for $n=2$, a $2\times 2$ matrix). Using Eqs. (2) and (3), one then obtains, by employing also its equilibrium counterpart $P^{\text{eq}}_{ij}(t)$,

The considered 4-state process is exactly solvable (see Supplemental Material Sup ), and $\mathcal{S}^{\prime}$, $\mathcal{D}^{\prime}$ and $P^{\text{eq}}$ so obtained are shown in Fig. 1(a). When applying the scheme experimentally, these curves are to be measured.

In this example, we take $O(X)=X$, i.e., we consider the response of $\langle X\rangle.$ The corresponding $\chi^{\text{rf}}$ are then found via Eq. (11) which, using $n=2$, simplifies to

Since $\mathcal{S}^{\prime}_{ij}$ is anti-symmetric and we have $n=2$, the sums reduce to the term $0\to 1$, and the nontrivial second order is the product of the functions shown in Fig. 1(a).

We show analytically Sup that Eq. (13b) indeed yields the exact second order response, which, having coarse grained a four state to a two state model, is an explicit confirmation of the proposed scheme.

Fig. 1 (b) shows the resulting $\chi_{2}$ as a function of time for an extreme choice of parameters: The rates $A\leftrightarrow B$ and $C\leftrightarrow D$ are small compared to the rates $B\leftrightarrow C$. Because of this, the density $\rho_{A}$ relaxes much slower than $\rho_{B}$ (inset), and the $\chi_{2}(t)$ shows two distinct time scales. This demostrates that Eq. (9) does neither rely on fast relaxation of integrated degrees, nor on Markovianity of the resulting two state system. For $t\to\infty$, $\chi_{2}$ vanishes because of symmetries.

2d Ising model – To demonstrate practical applicability, we consider Ising model on a periodic square lattice with nearest-neighbor interactions among N spins $s_{i}=\pm 1$ and following Metropolis dynamics Newman and Barkema (1999), studied using Monte-Carlo simulations ²²2A randomly selected spin flips with a rate $\min\{1,e^{-\Delta H/T}\}$ where $\Delta H$ is the change in the Hamiltonian due to the proposed flip.. See Ref. Lippiello et al. (2008b) for nonlinear response theory in the Ising model. The Hamiltonian

is asymmetric due to the presence of a magnetic field $h$ (included to allow for a finite $\chi_{2}$). $\varepsilon$ gives the strength of perturbation which acts on ${\cal N}\leq N$ spins, and the unitstep function $\Theta(t)=0$ if $t<0$ and $\Theta(t)=1$ otherwise. With $k_{B}=1$, $h$ and temperature $T$ are dimensionless. For $h=0$, the 2d Ising model shows a paramagnet-ferromagnet-transition at temperature $T_{c}\simeq 2.269$ Baxter (1982). Our finite system with a lattice of size $N=16\times 16$ and $T=2.45$ shows ferromagnetic order, however randomly flipping collectively the sign of the magnetization $m=\frac{1}{N}\sum_{i=1}^{N}s_{i}$, on a slow time scale.

For the macrovariable $X=\sum_{i=1}^{\cal N}s_{i}$, corresponding to $n={\cal N}+1$ macrostates, the perturbation in Eq. (14) is of the form $V(X)$ (namely $V(X)=X$). An extreme limit is a local perturbation (${\cal N}=1$), where only a single tagged spin is perturbed. Here, the interpolation scheme is applied by only tracking (measuring) the dynamics of that tagged spin ($n=2$), while the configuration of the surrounding spins need not be known ³³3We omit presentation of the numerical data for this case here..

More challenging, we consider a global perturbation (${\cal N}=N$), aiming at the sign of the magnetization as chosen observable of interest, specifically $O=\Theta(m)$. With $h=0.005$, $\langle O\rangle^{\rm eq}\simeq 0.613$ in the equilibrium state. Does one need $N+1=257$ macrostates in this case? Practically, a much smaller number turns out to be sufficient. We use $n=2$, $4$ and $6$ (see sketch in Fig. 2), ruling out odd values for symmetry.

In our simulations, we measure ${\cal S}^{\prime}$ and ${\cal D}^{\prime}$ with a small value of $\varepsilon=0.0005$ using Eq. (12) ⁴⁴4To increase precision, all curves are measured at $\pm\varepsilon$. The results for $\varepsilon$ and $-\varepsilon$ are added (linear in $\varepsilon$) and subtracted (second order in $\varepsilon$), making use of symmetry Sup ; Helden et al. (2016).. This yields the curves in Fig. 2(a) (for ease of presentation, we only show the case $n=2$). The predicted second order response, $\chi^{\text{rf}}_{2}(t)$, is then given by Eq. (11b), i.e., summing over the matrix elements of ${\cal S}^{\prime}$, ${\cal D}^{\prime}$ and $P^{\rm eq}$. For $n=2$, this sum in given in Eq. (13b), and contains only one term: It is the product of the functions in Fig. 2(a). For larger $n$, more terms are summed. This yields the curves in Fig. 2(b). We also measured the second order response using the conventional way (see Eq. (10b)), for which we have used a larger value of $\varepsilon=0.003;$ see the (blue) curve denoted ‘per’ in Fig. 2(b). The very good agreement in Fig. 2(b) confirms the main claim of the paper: We used simulations at $\varepsilon=0$ and $\varepsilon=0.0005$, and obtained the nontrivial extrapolation to a larger perturbation $\varepsilon=0.003$. As a practical aspect, the conventional way of determining $\chi_{2}$ (using Eq. (10b)) needs about ten times the amount of computational effort to obtain curves with similar statistics. The curves for different $n$ can only be distinguished in a logarithmic presentation (Fig. 2(b) inset), where the long time limit, found in a static measurement Sup , is indeed aproached better and better for increasing $n$. We note that for other systems, the convergence with $n$ may be slower.

The scheme amounts to measuring transitions rates between the different values of $X$ which are, in suitable systems, obtained much more easily compared to the measuments needed for microscopic response theory. Once experimental trajectories are obtained, the transition rates can be evaluated for different $n$, so that, larger $n$s do not necessarily require more experimental measuring time.

Recapitulating, $V=V(X)$ is a sufficient condition for accuracy of the proposed scheme. It means that unperturbed degrees of freedom can be coarse grained straightforwardly. In our examples, these are the unperturbed links or spins, but, in general, these can also include spatial or momentum degrees of freedom. Practically, we noted that the condition $V=V(X)$ is not necessary, so that much coarser descriptions as implied by this condition can suffice. By testing convergence with $n$, the accuracy of the method can be controlled. Because naturally, the obtained resolution for the observable is limited by the number of macrostates, this approach is especially useful if the behavior of a low dimensional observable is sought, such as the order parameter of a (phase) transition.

The presented coarse graining and extrapolation scheme constitute a conceptually new approach to nonlinear response theory. Because micro-degrees do not have to be monitored, it has a large range of applicability in complex systems. While circumventing the experimental need of applying strong perturbations, the scheme can also be more efficient regarding computation time compared to the traditional way of obtaining response functions, which is of additional advantage for slow systems. We note that at any order of perturbation, the response formula contains the change of $\mathcal{D}$ in at most one order lower, so that we expect the extrapolation to be extendable beyond second order.

Future work will investigate time-dependent perturbations, and perturbations via nonconservative force fields.

Supplemental Material for “Extrapolation to nonequilibrium from coarse grained response theory”

4-state jump process: The 4-state jump process [refer to Fig. 1 in the main text], provides an example for which the nonlinear response can be exactly calculated, and exact validity of the extrapolation scheme can be demonstrated. The coarse grained path probabilities $P_{ij}(t)$ are expressed as sum over the microscopic paths connecting the macrostates $i$ and $j.$ For example, the path with initial state $i=0$ and final state $j=1$ at $t$ has a probability,

where $p_{\alpha\beta}(t)$ denotes the probability that starting from micro-state $\alpha$ at $t=0$ the system reaches state $\beta$ in time $t$. We obtained it exactly by solving the time dependent Master equation for any choice of jump rates; $\rho_{\alpha}^{\text{eq}}$ is the equilibrium probability for the system to be in the micro-state $\alpha.$

We consider, as in Fig. 1 in the main text, a perturbation which changes the ‘forward’ jump rate connecting the two macrostates: $k_{BC}=e^{\varepsilon}.$ The matrices ${\cal S}^{\prime}$ and ${\cal D}^{\prime}$ are found form Eq. (8) in the main text, using $P_{ij}(t)$ from Eq. (LABEL:eq:4st_P). Taking ratios of the equilibrium and perturbed macro-path probabilities, we obtain explicitly,

where $\lambda=\sqrt{1+r^{2}}.$ The second order derivative, $\mathcal{S}^{\prime\prime}=0$, vanishes exactly. The second order response of $\langle X\rangle,$ is then [following Eq. (9) in the main text] predicted,

On the other hand, the second order response can also be found analytically without use of the response formula, by Taylor’s expansion of the exact nonequilibrium result $\langle X(t)\rangle=\rho_{C}(t)+\rho_{D}(t)$ where $\rho_{\alpha}(t)$ is the density in the micro-state $\alpha$ at time $t$ in the perturbed process. This yields

which matches exactly with Eq. (20) using (LABEL:eq:4stSD1), so that $\chi_{2}^{\text{eq}}=\chi_{2}$ is demonstrated.

Static response formula for second order response: In the case of a potential perturbation, the stationary long time values of the second order response $\chi_{2}$ can independently be obtained by a static response formula, which results from a Taylor’s expansion of the Boltzman weight of the perturbed system.

For the Ising model discussed in the main text [see Eq. (11) therein], and for the case ${\cal N}=N$ (i.e., the global perturbation discussed in the main text), second order response for observable $O$ is given by,

where $m$ is the magnetization; $\langle\cdot\rangle$ denotes expectation in the unperturbed state and $\langle A;B\rangle=\langle AB\rangle-\langle A\rangle\langle B\rangle$ denotes the connected correlation. The above expression is used to compute $\chi_{2}^{\text{st}}$ which is displayed in Fig. 2 of the main text. We observe that this method yields a rather precise result which is used as a long time benchmark for the time dependent solutions.

Improving accuray of response data: Fig. 2 in the main text compares different methods (different numbers of macrostates) for finding second order responses. In order to obtain accurate data for this comparison, and to minimize sources of errors beyond the coarse graining, we measured all perturbed states under perturbations of $\varepsilon$ as well as $-\varepsilon$. We denote the corresponding transition rates $P_{ij}^{\varepsilon}$ and $P_{ij}^{-\varepsilon}$. Eq. 9 in the main text transforms then to

These expressions have an error of order $\mathcal{O}(\varepsilon^{2})$, while Eq. 9 in the main text has an error $\mathcal{O}(\varepsilon)$. Eq. (24) has been used to evaluate $\mathcal{S^{\prime},D^{\prime}}$ to find the curves in Fig. 2 in the main text.

Similarly, the directly measured second order response for $\langle O\rangle$ has been obtained from

where $\langle O\rangle_{\varepsilon}$ and $\langle O\rangle_{-\varepsilon}$ denote the expectation value of $O$ at perturbation strengths $\varepsilon$ and $-\varepsilon,$ respectively. Also for Eq. (25), the error ($\mathcal{O}(\varepsilon^{2})$) is reduced compared to Eq. (7b) in the main text, which generaly yields an error of order $\mathcal{O}(\varepsilon)$.

Using Eqs. (24) and (25) thus improves accuray of the obtained curves in Fig. 2. We note that this method contains no principle change in strategy, and does not render the comparison of the proposed scheme to the conventional method, as it improves the data in both methods equally (changing errors from order $\mathcal{O}(\varepsilon)$ to $\mathcal{O}(\varepsilon^{2})$). The formulas presented in the main text yield good data as well, however, the difference between different numbers of macrostates $n$ would be less easily apparent. Discussing and testing the behavior of the extrapolation scheme as a function of $n$ is an important aspect of this manuscript, so that we seek data as accurate as possible.