Nonlinear transport coefficients from large deviation functions

TL;DR

This paper presents a novel method to compute high-order nonlinear transport coefficients in stochastic systems using large deviation theory, linking response to equilibrium correlations and enabling efficient evaluation.

Contribution

It introduces a thermodynamic-like relation for nonlinear response based on large deviation functions, applicable to complex nanoscale systems.

Findings

Successfully predicts transport coefficients in single particle systems.

Demonstrates efficiency in systems with thermal rectification.

Establishes a practical route for nonlinear response evaluation.

Abstract

Nonlinear response occurs naturally when a strong perturbation takes a system far from equilibrium. Despite of its omnipresence in nanoscale systems, it is difficult to predict in a general and efficient way. Here we introduce a way to compute arbitrarily high order transport coefficients of stochastic systems, using the framework of large deviation theory. Leveraging time reversibility in the microscopic dynamics, we relate nonlinear response to equilibrium multi-time correlation functions among both time reversal symmetric and asymmetric observables, which can be evaluated from derivatives of large deviation functions. This connection establishes a thermodynamic-like relation for nonequilibrium response and provides a practical route to its evaluation, as large deviation functions are amenable to importance sampling. We demonstrate the generality and efficiency of this method in predicting transport coefficients in single particle systems and an interacting system exhibiting thermal rectification.

Figures

Equations

$$\mathcal{J}=L_0+L_{1}X+L_2{X}^2+L_3{X}^3+\cdots ,$$

$$\frac{P_X[\tilde{x}]}{P_0[\tilde{x}]}=e^{\beta U[\tilde{x}]},$$

$$U[\tilde{x}]=A[\tilde{x}]+S[\tilde{x}],$$

$$A[\tilde{x}]=t_{N}JX,J[\tilde{x}]=\frac{1}{t_N}{\int }_0^{t_N}dtj({\tilde{x}}_t)$$

$$S[\tilde{x}]=\sum _{i=1}^p{t}_N{Q}_i{X}^i,Q_i[\tilde{x}]=\frac{1}{t_N}{\int }_0^{t_N}dt{q}_i({\tilde{x}}_t),$$

$${\psi }_X(\mathbf{\lambda })=\underset{t_N\to \infty }{\lim }\frac{1}{t_N}\ln {⟨e^{-({\lambda }_J{t}_{N}J+{\sum }_i{\lambda }_{Q_i}{t}_N{Q}_i)}⟩}_X,$$

$${\psi }_X({\lambda }_J,$$

$$={\psi }_0({\lambda }_J-\beta X,{\lambda }_{Q_1}-\beta X,{\lambda }_{Q_2}-\beta X^2,\cdots ),$$

$$\left\langle J\right\rangle_{X}=-\frac{\partial\psi_{X}({\bm{\lambda}})}{\partial\lambda_{J}}\Bigr{|}_{\begin{subarray}{c}\bm{\lambda}=\mathbf{0}\end{subarray}}=-\frac{\partial\psi_{0}({\bm{\lambda}})}{\partial\lambda_{J}}\Bigr{|}_{\begin{subarray}{c}\lambda_{J}=-\beta X,\lambda_{Q_{i}}=-\beta X^{i}\end{subarray}}.$$

$${⟨J⟩}_X$$

$$+\frac{{\beta }^2{t}_N^2}{2}[{⟨(\delta J{)}^3⟩}_0+{⟨\delta J(\delta Q_1{)}^2⟩}_0+2{⟨(\delta J{)}^2\delta Q_1⟩}_0$$

$$+\frac{2}{\beta t_N}{⟨\delta J\delta Q_2⟩}_0]X^2+\cdots ,$$

$${⟨J⟩}_X=\frac{{⟨J{e}^{\beta U}⟩}_0}{{⟨e^{\beta U}⟩}_0}.$$

$${⟨J⟩}_X=\beta t_N{⟨J^2⟩}_{0}X+{\beta }^2{t}_N^2{⟨J^2{Q}_1⟩}_0{X}^2+O(X^3).$$

$$L_{1}=\beta\frac{\partial^{2}\psi_{0}}{\partial\lambda^{2}_{J}}\Bigr{|}_{\begin{subarray}{c}\bm{\lambda}=\mathbf{0}\end{subarray}}\,,\quad L_{2}=-\beta^{2}\frac{\partial^{3}\psi_{0}}{\partial\lambda_{J}^{2}\partial\lambda_{Q_{1}}}\Bigr{|}_{\begin{subarray}{c}\bm{\lambda}=\mathbf{0}\end{subarray}}.$$

$$\gamma \dot{r}=-\frac{d{V}_p(r)}{dr}+F+\eta .$$

$$j({\tilde{x}}_t)=\dot{r}(t),q_1({\tilde{x}}_t)={\gamma }^{-1}\frac{d{V}_p}{dr}(t).$$

$${⟨J⟩}_F$$

$$+{\beta }^3{t}_N^3[\frac{1}{6}({⟨J^4⟩}_0-3{⟨J^2⟩}_0^2)+$$

$$\frac{1}{2}({⟨J^2{Q}_1^2⟩}_0-{⟨J^2⟩}_0{⟨Q_1^2⟩}_0)]\frac{F^3}{8}+O(F^4).$$

$$L_{\mathbf{\lambda }}=-\frac{\partial V_p}{\partial x}\left(\frac{\partial }{\partial x}-{\lambda }_J\right)+T{\left(\frac{\partial }{\partial x}-{\lambda }_J\right)}^2-{\lambda }_{Q_1}\frac{\partial V_p}{\partial x},$$

$$j({\tilde{x}}_t)=\frac{T_0}{4}\left(m_1{\dot{v}}_1{v}_1-v_1{F}_1^{\mathrm{FPUT}}-m_N{\dot{v}}_N{v}_N+v_N{F}_N^{\mathrm{FPUT}}\right),$$

$$P_{\Delta T}[\tilde{x}]$$

$$\times \prod _l\sqrt{\frac{m_N}{2\pi T_N}}\exp [-\frac{m_N{v}_N(t_l{)}^2}{2T_N}]$$

$$Q_1[\tilde{x}]=\frac{T_0}{4t_N}\left(\sum _k{m}_1{v}_1^2(t_k)-\sum _l{m}_N{v}_N^2(t_l)\right)-J[\tilde{x}].$$

$${⟨J⟩}_{\Delta T}$$

$$+{\beta }^2{t}_N^2{⟨J^2{Q}_1⟩}_0{\left(\frac{2\Delta T}{T_0^2}\right)}^2+O(\Delta T^3)$$

$${\dot{v}}_x=-\gamma (v_x-U_x)+F+{\eta }_x,{\dot{v}}_y=-\gamma (v_y-U_y)+{\eta }_y,$$

$$j({\tilde{x}}_t)=v_x(t),q_1({\tilde{x}}_t)={\dot{v}}_x(t)/\gamma -U_x(t).$$

$${⟨J⟩}_F=\beta \left(\frac{1}{\gamma \beta }+\frac{t_N}{2}{⟨J^2⟩}_0-\frac{t_N}{2}{⟨Q_1^2⟩}_0\right)\frac{F}{2}+O(F^2).$$

$$m_i{\dot{v}}_i=-m_i\gamma v_i-\frac{dV(\mathbf{r})}{d{r}_i}+F+{\eta }_i,$$

$$P_F[\tilde{x}]\propto$$

$$\exp \left[-\int dt\sum _{i=1}^N\frac{(m_i{\dot{v}}_i+m_i\gamma v_i+dV(\mathbf{r})/d{r}_i-F{)}^2}{4m_i\gamma k_{B}T}\right].$$

$$U[\tilde{x}]=-\sum _{i=1}^N\frac{-(2m_i{\dot{v}}_i+2m_i\gamma v_i+2dV(\mathbf{r})/d{r}_i)F+F^2}{4m_i\gamma }.$$

Full text

irreversible process $|$ nanoscale transport $|$ nonequilibrium

††preprint: 1

Transport processes are often described by phenomenological laws, such as Fourier’s law for heat conduction, Fick’s law for diffusion, and Ohm’s law for electrical conduction. These equations propose a relationship between the thermodynamic force, or affinity, that drives the system out of equilibrium and the resultant current, where the strength of the response is described by a material specific transport coefficient. While connections between these coefficients and the underlying microscopic degrees of freedom are well established when forces are small, such connections are much less clear when forces are large and responses are nonlinear. The purpose of this paper is to develop a general scheme for the calculation of arbitrarily high order transport coefficients in stochastic systems from microscopic equilibrium fluctuations, and illustrate its viability by computing nonlinear transport coefficients from the derivatives of large deviation functions.

The notion of relating transport coefficients to spontaneous fluctuations in equilibrium dates back to Onsager’s regression hypothesis, which states that the relaxation of a system from a given nonequilibrium state prepared by an external perturbation, obeys the same law as the regression of the same state produced by spontaneous thermal fluctuations Onsager (1931). A further connection between the two was established later by deriving the path probability of a given succession of states in a system where the macroscopic variables are Gaussian distributed Onsager and Machlup (1953). Along the same line, Kubo constructed a general statistical-mechanical theory for irreversible processes Kubo (1957). The well known Green-Kubo formulas express a linear response coefficient with the integral of the current autocorrelation function over time, and have become a standard method to compute linear transport coefficients Green (1954); Levesque et al. (1973); Schelling et al. (2002). Onsager and Kubo’s work are the cornerstones of linear response theory, which applies to irreversible processes close to equilibrium.

As the external perturbation grows stronger, nonlinear transport behaviors naturally arise, such as shear thinning in complex fluids, and current rectification in electrical transistors. Such behaviors become especially prevalent in nanoscale systems, where microscopic fluctuations play an important role. As a consequence, understanding the origin of nonlinear effects are key to the design and manipulation of nanoscale devices. One of the most efficient methods to compute nonlinear transport coefficient is direct nonequilibrium molecular simulation, where a current is driven through the system by application of specific boundary conditions Tenenbaum et al. (1982); Baranyai and Cummings (1999), or by altering the equations of motion Hoover et al. (1980); Evans (1982); Müller-Plathe (1997). Such methods, however, are generally not transferable among different transport processes, and the result can be sensitive to how the current is generated Tuckerman et al. (1997); Tenney and Maginn (2010). Moreover, with direct simulation the computed response is not easily connected to specific molecular degrees of freedom.

Kubo in his original framework also derived a perturbation expansion for the nonequilibrium phase space distribution. However, such expression, as well as many that followed Efremov (1969); Stratonovich (1970), have been difficult to translate into measurable forms with clear physical interpretations. Several non-perturbative treatments were proposed around the same time Yamada and Kawasaki (1967); Bochkov and Kuzovlev (1977), yet an explicit expression for computing nonequilibrium averages is lacking. The simplest generalization of the Green-Kubo formulas that fills in this gap is perhaps the transient time correlation function formalism Visscher (1974); Dufty and Lindenfeld (1979); Cohen (1983). It gives an exact relation between the nonequilibrium ensemble average of an observable and the transient time correlation function between the observable and the current. In practice, a nonequilibrium simulation has to be performed, yet the statistical uncertainty in the nonequilibrium average is often much larger than that of the direct nonequilibrium simulation Evans and Morriss (1988).

The last two decades have observed a growing interest in the field of stochastic thermodynamics, stimulated by new possibilities in experimental interrogation of nanoscale systems. Nonequilibrium fluctuation theorems Jarzynski (1997); Crooks (1999); Gallavotti and Cohen (1995) and thermodynamic uncertainty relations Barato and Seifert (2015); Gingrich et al. (2016) impose constraints on the distribution functions of fluctuating thermodynamic quantities like heat and work, which represent refinements of the second law to driven systems. These theoretical developments have enabled the establishment of a more formal theory for nonequilibrium steady states. In particular, a considerable amount of work has been devoted to predicting response coefficients for stochastic systems. Extended fluctuation-dissipation theorems have been derived for linear responses of nonequilibrium steady states Harada and Sasa (2005); Speck and Seifert (2006, 2009); Baiesi et al. (2009), which can be translated into second order response around equilibrium Basu and Maes (2015). More generally, multivariate fluctuation relations imply connections between transport coefficients and cumulants of the current Andrieux and Gaspard (2004); Gaspard (2013). Many of these advances have been enabled by large deviation theory Ellis (2007); Touchette (2009), which provides a set of mathematical tools for characterizing and evaluating fluctuations in nonequilibrium systems. These tools have been employed widely to illustrative model systems, which not only serve as testing grounds for theories, but also reveal new insights into nonequilibrium phenomena Hedges et al. (2009); Hurtado and Garrido (2011); Limmer and Chandler (2014); Pre and Limmer (2018).

In the present work, we aim to translate these earlier efforts to a larger class of nonequilibrium systems, and derive an expression for arbitrarily high order transport coefficients in terms of equilibrium correlation functions. To compute the correlation functions, we rely on the direct evaluation of large deviation functions (LDFs). We have shown in a previous work Gao and Limmer (2017) that linear transport coefficients can be extracted from the second derivative of LDFs in a statistically efficient way. Here, we extend the methodology to nonlinear transport regimes. In what follows, we first formulate our theory and then demonstrate our method in a collection of examples, including a Brownian ratchet, a model thermal rectifier, and a tracer particle in a convective flow.

Theory and Methodology

We are interested in stochastic systems maintained in nonequilibrium steady states by an affinity $X$, which can either be an external force, or contacts with different thermodynamic reservoirs. We assume that the macroscopic current $\mathcal{J}$ generated in the system as a response to $X$ can be written as a polynomial expansion,

[TABLE]

where $L_{0}$ is the current in the absence of the affinity, which should vanish for equilibrium systems. $L_{1}$ is known as the linear transport coefficient, and all the higher order coefficients $L_{2},L_{3},\cdots$, are nonlinear transport coefficients. The generality of this polynomial expansion will be discussed towards the end of this paper. Given this assumption, in this section we derive an explicit expression for these transport coefficients. Although we will restrict ourselves to a single type of affinity, the generalization to multiple ones is straightforward.

Derivation of nonlinear coefficients from LDFs

We consider a continuous stochastic trajectory of length $t_{N}$ denoted by $\tilde{x}$, where $\tilde{x}_{t}=\{\bm{r}(t),\bm{v}(t)\}$ is the specific configuration of the system at time $t$, with coordinates $\bm{r}(t)$ and velocities $\bm{v}(t)$. We define a relative stochastic action $\beta U[\tilde{x}]$ by a ratio of the probability of observing a specific path with and without the affinity,

[TABLE]

where the subscripts denote the value of the affinity; $\beta=1/k_{B}T$ where $T$ is the temperature, $k_{B}$ is Boltzmann constant that will be set to 1 in the following calculations.

We next decompose the total relative path action into two parts according to their time reversal symmetry,

[TABLE]

with $A[\tilde{x}]$ denoting the asymmetric part, and $S[\tilde{x}]$ the symmetric part. Specifically, if we define an operator $\mathbb{T}$ which returns the time-reversed counterpart of a path $\tilde{x}$ so that $\mathbb{T}\tilde{x}_{t}=\{\bm{r}(t_{N}-t),-\bm{v}(t_{N}-t)\}$, it follows that $A[\mathbb{T}\tilde{x}]=-A[\tilde{x}]$, $S[\mathbb{T}\tilde{x}]=S[\tilde{x}]$. The time asymmetric part is the entropy production of the irreversible process, which can be written as the product of the affinity and the conjugated time extensive current Onsager (1931),

[TABLE]

The symmetric part is often referred to as the activity Maes et al. (2008),

[TABLE]

The number of terms $p$ depends on system details, though for Gaussian processes it is often only 1. For details, we show the derivation of $J$ and $Q_{i}$’s for a general underdamped Langevin system in the Appendix.

Next, we construct a LDF of the form

[TABLE]

where $\bm{\lambda}=(\lambda_{J},\lambda_{Q_{1}},\lambda_{Q_{2}},\cdots)$, i.e. each dynamical variable is exponentially biased by a conjugated $\lambda$. The bracket denotes path ensemble average with respect to the path probability density $P_{X}[\tilde{x}]$. It follows by definition that

[TABLE]

which provides a symmetry that links the statistical bias ${\bm{\lambda}}$ to the physical driving $X$. This relation is distinct from a fluctuation theorem and encodes the fact that $X$ is not a conjugate variable to $J$, but rather a linear combination of $J$ and the $Q_{i}$’s. The derivatives of the cumulant generating function provide information about the self and cross correlation functions of $J$ and the $Q_{i}$’s. For example, the nonequilibrium average current is given by the first derivative

[TABLE]

Expanding $\psi_{0}({\bm{\lambda}})$ assuming $\beta X$ is small,

[TABLE]

we find a microscopic relation between the average current and affinity, where for brevity, we have only written down explicitly terms up to $O(X^{2})$. The notation $\left\langle A_{1}A_{2}\cdots A_{n}\right\rangle=t_{N}^{-n}\int_{0}^{t_{N}}d{\bf t}^{n}\left\langle a_{1}(\tilde{x}_{t_{1}})a_{2}(\tilde{x}_{t_{2}})\cdots a_{n}(\tilde{x}_{t_{n}})\right\rangle$ and $\delta A=A-\left\langle A\right\rangle_{0}$ are adopted throughout. By comparing with Eq. 1, and assuming the time averaged current is equal to the macroscopic current, the expansion in Eq. 9 provides explicit expressions for arbitrarily high order transport coefficients in terms of multi-time correlation functions in the absence of the affinity, which can be computed from $\psi_{0}({\bm{\lambda}})$. While multi-time correlation functions are in general difficult to converge by direct evaluation, $\psi_{0}({\bm{\lambda}})$ can be computed efficiently using importance sampling methods Nemoto and Sasa (2014); Nemoto et al. (2017); Klymko et al. (2018); Ray et al. (2018).

To highlight the novelty of our method, we discuss how our expression is different from a few previous results in the literature. Firstly, our result is intimately related to work by Maes et al Colangeli et al. (2011); Basu et al. (2015), though our current expression Eq. 9 is different beyond the third order. This stems from the fact that the current is derived in terms of a cumulant expansion instead of a moment expansion, as seen clearly in Eq. 9, and from rewriting Eq. 8,

[TABLE]

The cumulant expansion ensures better convergence for the transport coefficients, especially for dynamical variables whose distributions are close to Gaussian Freed (1968); Wang et al. (1998). Secondly, our method is distinct from the multivariate fluctuation relations that have been derived before by constructing the LDF of only current observables Gaspard (2013), where transport coefficients are expressed as mixed derivatives of both $\lambda$’s and the affinities. Here, by introducing the time symmetric observables $Q_{i}$’s, the knowledge of the equilibrium function $\psi_{0}({\bm{\lambda}})$ completely determines the current-affinity relationship.

Simplification by symmetries

Decomposing the total action by time reversal symmetry enables us to greatly simplify Eq. 9. If the reference system without the affinity is in equilibrium, then it obeys microscopic time reversibility $P_{0}[\mathbb{T}\tilde{x}]=P_{0}[\tilde{x}]$. It follows that the time reversal odd terms, such as the average current $\left\langle J\right\rangle_{0}$, and terms like $\left\langle JQ_{i}\right\rangle_{0}$, $\left\langle JQ_{i}^{2}\right\rangle_{0}$ and $\left\langle J^{3}\right\rangle_{0}$, will vanish. This reduces Eq. 9 to

[TABLE]

Note that the linear response reduces to the normal fluctuation-dissipation theorem, while higher order responses are described by higher order correlations between $J$ and $Q_{i}$’s. These correlation functions are given by derivatives of $\psi_{0}(\bm{\lambda})$. Specifically, the linear and first nonlinear transport coefficients are

[TABLE]

Further simplifications can be made if the system exhibits certain spatial symmetries so that the dynamics is unchanged upon inverting the coordinates along a specific axis.

Numerical evaluation of large deviation functions

To evaluate the large deviation function defined in Eq. 6, we use a diffusion Monte Carlo method known as the cloning algorithm Giardina et al. (2006), where an ensemble of trajectories is propagated in parallel. Each individual trajectory is known as a walker, and collectively they undergo a population dynamics that samples a biased trajectory ensemble. A more detailed description of the algorithm can be found in our earlier work Gao and Limmer (2017).

One of the key factors in determining the statistical efficiency of the cloning algorithm is the number of independent trajectories sampled over $t_{N}$. In all the calculations shown in this paper, walker numbers are chosen individually for each parameter so that at the end of the simulation, the number of independent walkers which have not been replaced is at least on the order of $10^{2}$.

While transport coefficients can also be computed directly through nonequilibium simulation methods, we note that by evaluating them from equilibrium fluctuations, our method is less sensitive to finite size effects arising from boundary conditions and altered equations of motion required to simulate a driven system. In addition, while direct simulation measures the nonequilibrium response at a finite value of the affinity, our method generically generates the response for a continuum of affinities. Indeed, if one is specifically interested in the current response at a finite field $X$, it suffices to measure the local curvature of $\psi_{X}$ at $\bm{\lambda}=\bm{0}$ by computing the equilibrium large deviation function $\psi_{0}$ around $\lambda_{J}=\beta X,\lambda_{Q_{i}}=\beta X^{i}$ according to Eq. 7, where we have established a direct connection between the statistical biased and nonequilibrium ensembles.

Models

Brownian Ratchet

To validate our method, we first consider a single Brownian particle moving on a one-dimensional asymmetric potential landscape $V_{p}(r)=\sin(r/l+\sin(r/l)/2)$. This is a simplest continuous model that exhibits an asymmetric response of the particle’s displacement to an additional constant force. The dynamics of the particle with unit mass are described by the overdamped Langevin equation

[TABLE]

where $F$ is a constant force, $\gamma$ is the friction and $\eta$ is the white noise satisfying $\left\langle\eta(t)\right\rangle=0,\left\langle\eta(0)\eta(t)\right\rangle=2\gamma\delta(t)/\beta$. We set $\gamma$ and $l$ to 1, to define a reduced unit system and take $\beta=1/0.3$. Considering the mobility of the particle in response to the affinity $X=F/2$, the dynamical variables take the form

[TABLE]

Here, the time symmetric part is characterized by the force exerted onto the particle by the ratchet potential. There is an additional second order term $q_{2}(\tilde{x}_{t})$, which is a constant and does not enter into the expression for the current. Expanding Eq. 9 up to the third order, we arrive at an expression for the integrated nonequilibrium current, or displacement,

[TABLE]

It is worth noting that the higher order responses are described not only by the current fluctuations, but also the correlations between current and force fluctuations. If the potential is symmetric, the second order term will vanish due to inversion symmetry. It is the asymmetry in the ratchet potential, and as a consequence, the asymmetry in the forces, that gives rise to the even order terms. Specifically, the correlation between the squared current and the force dictates the size of the rectification of the ratchet.

In this simple case, the LDF $\psi_{0}({\bm{\lambda}})$ can be computed numerically exactly by diagonalizing the tilted generator Touchette (2009),

[TABLE]

the largest eigenvalue of which is the LDF. To solve for the eigenvalues, we construct the tilted operator with a normalized Fourier basis set $\exp({ikx})$ where $k\in[-15,15]$ and is an integer. The 2D surface of $\psi_{0}({\bm{\lambda}})$ is shown in Fig. 1$A$. The obvious deviation from the normal distribution, especially the asymmetry in the $\lambda_{Q_{1}}$ direction, leads to the observed nonlinear behavior of $\langle J\rangle_{F}$. The correlation functions appearing in the average current are computed from numerical derivatives of $\psi_{0}$ as in Eq. 12, yielding $t_{N}\left\langle J^{2}\right\rangle_{0}=0.0184$, $t_{N}^{2}\left\langle J^{2}Q_{1}\right\rangle_{0}=-0.0194$, $t_{N}^{3}(\left\langle J^{4}\right\rangle_{0}-3\left\langle J^{2}\right\rangle_{0}^{2})=0.712$, and $t_{N}^{3}(\left\langle J^{2}Q_{1}^{2}\right\rangle_{0}-\left\langle J^{2}\right\rangle_{0}\left\langle Q_{1}^{2}\right\rangle_{0})=-0.033$. Fitting errors in the coefficients are negligible. The nonequilibrium simulation results are obtained by integrating the equation of motion using a second-order Runge-Kutta algorithm Honeycutt (1992) with a time step $h=10^{-3}$. Numerical results shown are averaged over $10^{3}$ realizations with a total observation time $t_{N}=10^{3}$. In Fig. 1$B$, we compare direct nonequilibrium simulation results with our theoretical predictions, which improve as we include higher order terms. It is worth noting that in this case, as the large deviation function can be solved conveniently by matrix diagonalization, the computational efficiency of our method is far superior compared to direct simulation, which requires sampling of the nonequilibrium dynamics.

Model Thermal Rectifier

We next turn to an interacting system, a 1D thermal rectifier. A thermal rectifier is a type of material with intrinsic structural asymmetry so that it exhibits an asymmetric heat transport response when a temperature gradient is applied. Here we model a 1D thermal rectifier using the linearly mass-graded Fermi-Pasta-Ulam-Tsingou (FPUT) chain, which has been shown to capture thermal rectification behavior Yang et al. (2007). The mass of the $i$th particle is $m_{i}=m_{\mathrm{min}}+(i-1)(m_{\mathrm{max}}-m_{\mathrm{min}})/(N-1),\hskip 5.69054pti=1,\cdots,N$, with $m_{\mathrm{max}}=20m_{\mathrm{min}}$, and $N=20$ is the total number of particles. The particles oscillate around their average position $b_{i}=ia$, while interacting with neighboring particles through the quartic FPUT potential, $V_{\mathrm{FPUT}}(\bm{r})=\sum_{i=1}^{N+1}\kappa(r_{i}-r_{i-1}-a)^{2}/2+c(r_{i}-r_{i-1}-a)^{4}/4,$ where $\kappa$ and $c$ are the harmonic and anharmonic coupling, and $a$ is the lattice constant. Fixed boundary conditions are applied by fixing two fictitious particles at $b_{0}$ and $b_{N+1}$. We set $m_{\mathrm{min}}=\kappa=a/2=1$, to define a dimensionless unit system, in which we also let $c=1$. The equation of motion is integrated by the velocity Verlet algorithm with a timestep of $h=5\times 10^{-3}$.

To apply an external temperature gradient, we place the particles on the two ends in contact with infinitely large thermal reservoirs kept at temperatures $T_{1}=T_{0}+\Delta T$ and $T_{N}=T_{0}-\Delta T$, respectively, where the average temperature $T_{0}=0.1$. This is implemented by an Andersen thermostat Frenkel and Smit (2001) on each of the two end particles, where the time interval $\Delta t$ between successive collisions are distributed as $P(\Delta t)=\Gamma e^{-\Gamma\Delta t}$ and $\Gamma=0.8$ is the coupling strength to the bath. The affinity for the heat transfer process is $X=1/(T_{0}-\Delta T)-1/(T_{0}+\Delta T)\approx 2\Delta T/T_{0}^{2}$. In the steady state limit, the heat current from the two ends must be same in magnitude but opposite in direction, so we define the heat current as

[TABLE]

where $F_{i}^{\mathrm{FPUT}}=-\partial V_{\mathrm{FPUT}}(\bm{r})/\partial r_{i}$. This can be interpreted as the rate of change in the kinetic and potential energy, averaged over contributions from particle $1$ and $N$, scaled by a factor of $T_{0}/2$. This constant factor is essential in matching our definition of the current.

For the Andersen thermostat, the path probability can be written as

[TABLE]

where $t_{k}$($t_{l}$) are times at which the first(last) particle collides with the bath. From this we can derive the relative path action, and the time-symmetric part follows as

[TABLE]

Note that since we have chosen to express the action in terms a Taylor expansion of $\Delta T$, we have an infinite series of time-symmetric parts; however, we could have formulated in terms of $\Delta(1/T)$, in which case there will only be one single term $Q_{1}$. The average current as a function of $\Delta T$ becomes

[TABLE]

The first term is the standard Green-Kubo result for the thermal conductivity in terms of an integrated heat flux autocorrelation function. The second term correlates the squared heat flux with the instantaneous temperature difference, and results in an asymmetric response of the current. Since the masses of the two ends are different, the time dependent temperature fluctuations on either side of the chain need not be the same. This expression illustrates that thermal current rectification is a product of microscopic correlations between instantaneous temperature gradients and heat fluxes.

The evaluation of $\psi_{0}({\bm{\lambda}})$ from diffusion Monte Carlo is shown in Fig. 2$A$, while the nonequilibrium simulation results of the current is shown in Fig. 2$B$. Both the LDFs and the nonequilibrium simulation results averaged over 9600 realizations are calculated for trajectories with $t_{N}=2\times 10^{5}$. The LDFs are evaluated at $\lambda_{J}\in[-32,32]$ for $\lambda_{Q_{1}}=10$, and $\lambda_{J}\in[-20,20]$ for $\lambda_{Q_{1}}=0,-10$. Eight independent samples are calculated at each combination of $\bm{\lambda}$, and standard deviations are plotted as error bars in Fig. 2$A$. For each set of samples at a specific $\lambda_{Q_{1}}$, a parabola is fit, and the error bars in inset of Fig. 2$A$ are the standard deviation among the 8 curvatures. To estimate the statistical error in the second order transport coefficient, the fitted slope in the inset of Fig. 2$A$ is evaluated individually for each of the 8 sample sets, and standard error of the mean is reported in Fig. 2$B$. While $t_{N}\left\langle J^{2}\right\rangle_{0}$ is measured by the curvature of $\psi_{0}$ with $\lambda_{Q_{1}}=0$, $t_{N}^{2}\left\langle J^{2}Q_{1}\right\rangle_{0}$ is measured by how the curvature changes as we change $\lambda_{Q_{1}}$. Given the definition of $Q_{1}$, the change of the curvature with $\lambda_{Q_{1}}$ directly reports the change in the thermal conductivity with a temperature gradient. To make a better comparison, we plot the shifted LDFs $\psi_{0}^{*}(\lambda_{J},\lambda_{Q_{1}})=\psi_{0}(\lambda_{J},\lambda_{Q_{1}})-\psi_{0}(\lambda_{J},\lambda_{Q_{1}}=0)$ so that all the curves have the same minimum value at $\lambda_{J}=0$. Our method correctly predicts the rectification behavior with high accuracy.

To demonstrate the statistical efficiency of our approach, in Fig. 3 we compare the statistical error in the evaluation of $t_{N}^{2}\left\langle J^{2}Q_{1}\right\rangle_{0}$ by direct evaluation and using the cloning algorithm. The direct evaluation, plotted as $\lambda_{J}=0$ in Fig. 3, is evaluated by computing the triple correlation function by brute force averaged among $1.2\times 10^{5}$ independent trajectories. In the cloning algorithm, we evaluate $\psi(\lambda_{J},\lambda_{Q_{1}}=10)$ for various $\lambda_{J}$’s with the same number of walkers $N_{w}=1.2\times 10^{5}$, and compute the correlation function by $t_{N}^{2}\left\langle J^{2}Q_{1}\right\rangle_{0}=2[\psi_{0}(\lambda_{J},\lambda_{Q_{1}})-\psi_{0}(\lambda_{J},0)-\psi_{0}(0,\lambda_{Q_{1}})]/(\lambda_{J}^{2}\lambda_{Q_{1}})$. Statistical errors are estimated by the standard deviation from 9 independent simulations. We have chosen the same number of walkers and independent trajectories to ensure that the computational effort in terms of number of integration steps is the same. For equal computational effort evaluating $\left\langle J^{2}Q_{1}\right\rangle_{0}$ from the large deviation function exhibits smaller statistical error for sufficiently large $\lambda_{J}$’s, as shown in Fig. 3. For a fixed $\lambda_{Q_{1}}$, the error in the curvature of $\psi_{0}$ with respect to $\lambda_{J}$ should scale as $\propto\sqrt{\psi_{0}^{\prime\prime}(\lambda_{J})/\tilde{N}_{w}}/\psi_{0}(\lambda_{J})=1/(\lambda_{J}^{2}\sqrt{\tilde{N}_{w}})$, where $\tilde{N}_{w}$ is the number of uncorrelated walkers. Even though the correlation between the walkers as $\lambda_{J}$ increases brings in a non-trivial dependence of $\tilde{N}_{w}$ on $\lambda_{J}$, overall the cloning algorithm still out-performs the direct evaluation for $|\lambda_{J}|>8$ by about an order of magnitude.

Tracer particle in a convective flow

In the last section, we discuss a system out of equilibrium even at $X=0$. This is fundamentally different from the two cases above, as even at linear response traditional Green-Kubo formulas do not apply. Nevertheless the response of the system is still encoded in a LDF computed at $X=0$, though one evaluated in the nonequilibrium steady state. The previously derived extended fluctuation-dissipation theorems Harada and Sasa (2005); Speck and Seifert (2006, 2009) follow naturally from our general expression for the nonequilibrium current in Eq. 9.

We consider an underdamped particle with unit mass moving in a 2D hydrodynamic flow

[TABLE]

where $U_{x}=\partial\phi(\bm{r})/\partial r_{y}$ and $U_{y}=-\partial\phi(\bm{r})/\partial r_{x}$ describe the divergenceless flow of the stream function $\phi(\bm{r})=LU_{0}\sin(2\pi r_{x}/L)\sin(2\pi r_{y}/L)/(2\pi)$. We set $U_{0}=L=1$, which sets a scale for length and time. The noise $\eta_{i}$’s satisfy $\left\langle\eta_{i}(t)\right\rangle=0$, $\left\langle\eta_{i}(0)\eta_{j}(t)\right\rangle=2\gamma\delta_{ij}\delta(t)/\beta,i=x,y$. When the constant force $F=0$, it is the non-gradient form of the velocity field that forces the system out of equilibrium, though the average particle current still vanishes due to spatial symmetry. As a consequence, the linear response of the particle current includes two terms proportional to $X$, as shown in Eq.  9. The extended fluctuation-dissipation theorem includes a term that is the correlation between the current and a time symmetric variable derived from the first order part of the relative action as has been described previously Baiesi et al. (2009). Defining the affinity $X=F/2$, the dynamical variables are

[TABLE]

The time symmetric term includes the particle’s inertia, relative to its local flow velocity, which is only a function of the particle’s position. The linear response of the current can be viewed as a Green-Kubo relation, but correlating particle’s velocity relative to the velocity field $U_{x}$ Harada and Sasa (2005); Speck and Seifert (2006), up to a boundary term from the integral of $\dot{v}_{x}$.

We note that the dynamical variables satisfy the relation $j(\tilde{x}_{t})+q_{1}(\tilde{x}_{t})=\eta_{x}/\gamma$ when $F=0$, which implies that $t_{N}^{n-1}\left\langle(J+Q_{1})^{n}\right\rangle_{0}=(2/(\gamma\beta))^{n/2}(n-1)!!$ if $n$ is even, and vanishes if $n$ is odd. Note that this type of relationship is not restricted to this specific model - it is quite general for stochastic systems with quadratic path actions. These constraints on the moments reduce the number of unknown moments we have to compute. Here, we use the fact that the second order moments satisfy $t_{N}\left\langle(J+Q_{1})^{2}\right\rangle_{0}=2/(\gamma\beta)$, which allows us to rewrite the linear response as

[TABLE]

To compute the linear response coefficient, all we need are the curvatures of $\psi_{0}({\bm{\lambda}})$ along $\lambda_{J}=0$ and $\lambda_{Q_{1}}=0$, which are easily computable from diffusion Monte Carlo Giardina et al. (2006); Gao and Limmer (2017).

In the following calculations, we set $\gamma=0.1$, $\beta=0.5\times 10^{4}$. The underdamped equation is integrated with a second order Verlet-like integrator Sivak et al. (2013), with a timestep of $h=10^{-3}$. Both nonequilibrium results and LDFs are calculated from trajectories of the length $t_{N}=8\times 10^{5}$. Fig. 4$A$ and $B$ shows the LDF while fixing $\lambda_{Q_{1}}=0$ and $\lambda_{J}=0$, respectively. To evaluate $t_{N}\left\langle J^{2}\right\rangle_{0}$, we compute the LDFs at $\lambda_{J}\in[-0.0035,0.0035]$, $\lambda_{Q_{1}}=0$ and fit the curve with a parabola to estimate its curvature. Similarly, $t_{N}\left\langle Q_{1}^{2}\right\rangle_{0}$ is estimated with LDFs at $\lambda_{J}=0$, $\lambda_{Q_{1}}\in[-0.0035,0.0035]$. The statistical errors in LDFs are estimated by the standard deviation among 15 independent samples. To evaluate the statistical error in the linear transport coefficient, we fit a parabola to each of the sample sets, and compute the standard error of the mean among the 15 curvatures. The fitted curves yield an estimate of $t_{N}\left\langle J^{2}\right\rangle_{0}/2=2.6874\pm 0.0015$, $t_{N}\left\langle Q_{1}^{2}\right\rangle_{0}/2=2.6773\pm 0.0013$, and the mobility is $30.2\pm 6.8$.

In Fig. 4$C$, we plot the theoretically predicted linear response along with the nonequilibrium simulation results to show their agreement at small values of $F$. Nonequilibrium results are averaged over $2400$ independent trajectories and standard errors of the mean are plotted. This model has recently been shown to exhibit a negative differential mobility Sarracino et al. (2016). Near equilibrium, the mobility is proportional to $\langle J^{2}\rangle_{0}$, and thus must be non-negative. However, the linear response around a nonequilibrium steady state given in Eq. 23 clarifies how a negative differential mobility is possible. While not true at the conditions we consider, in principle $t_{N}\langle Q_{1}^{2}\rangle_{0}$ may be larger than $2/(\gamma\beta)+t_{N}\left\langle J^{2}\right\rangle_{0}$, resulting in a current that decreases with added force.

Discussion and Conclusion

We have introduced a method to compute arbitrarily high order transport coefficients that is heavily built on early works by Onsager and Kubo, and can be seen as an extension to stochastic systems in the modern language of stochastic thermodynamics. The relative action we introduce acts as the generalized “dissipation function” as Onsager first invoked Onsager and Machlup (1953), which gives the probability of a temporal succession of states. While in the linear response regime, the dissipation function is simply given by the rate of entropy production, to describe nonlinear response, additional information, namely the symmetric part of the action, is required. By constructing a LDF of both the time symmetric and asymmetric dynamical variables, we arrive at a function evaluated in equilibrium that contains all of the information about current-affinity relationship as encoded through microscopic correlations.

A primary limitation of our method is the assumption that the nonequilibrium phase space distribution function can be expanded as a power series of the affinity. Critics of Kubo’s nonlinear response theory have pointed out that for many transport processes, such expansions do not exist Evans (1983); Morriss and Evans (2013). Indeed, a perturbative treatment breaks down when long-range or long-time correlations exist in a system in the thermodynamic limit. This can happen in the vicinity of phase transitions, where diverging correlation lengths will cause a divergence in the correlation functions corresponding to the transport coefficients. As a result, our simple series expansion is likely to fail in describing dynamical phase transitions between nonequilibrium steady states. Additionally, constraints on dynamics can also result in long-range correlations. For example, low dimensional molecular fluids are known to exhibit diverging diffusivity computed from the Green-Kubo formula, which originates from the slowly decaying correlation of hydrodynamic modes associated with conserved quantities Alder and Wainwright (1970); Dorfman et al. (1994). In deterministic systems, conservation of energy and momentum confines the trajectories of the system to certain manifolds of the phase space. As a consequence, currents associated with conserved quantities are often non-analytic functions of the affinity Kawasaki and Gunton (1973); Yamada and Kawasaki (1975). To circumvent such problems, we have restricted ourselves to Markovian stochastic systems in the present work, where ergodicity is guaranteed. In all the examples shown, correlation functions decay to zero in microscopic timescales, so that their integrals in the long time limit are well defined.

The method we propose is general enough that it can be applied to stochastic systems both near and far from equilibrium, regardless of the specific type of affinity, as long as a microscopic definition of current and path action can be written down. Unlike direct nonequilibrium simulation, our results do not depend on the details of the system, such as the nature of the noise, since correlation functions are evaluated in an equilibrium ensemble. Our method is expected to exhibit superior statistical performance compared to methods relying on transient time correlation functions, based on a detailed study on the statistical error in the computation of linear transport coefficients by the Green-Kubo formulas and our method Gao and Limmer (2017). Higher order transport coefficients demand a more accurate evaluation of the large deviation function. However, with the advancement of algorithms enhanced by importance sampling techniques, we expect that our method will become a standard approach to computing higher order transport coefficients. Furthermore, the expression of transport coefficients in terms of multi-time correlation functions builds a connection between macroscopic transport processes and microscopic observables. We expect that theoretical manipulations on these multi-time correlation functions will give us further insights into how nonlinear response behaviors arise from the molecular details of the system.

Appendix

Derivation of the nonequilibrium average current for a general Langevin system

Consider a system with $N$ degrees of freedom interacting through potential $V(\bm{r})$ at temperature $T$. Each degree of freedom $i$ evolves under the underdamped Langevin equation

[TABLE]

where $m_{i}$ is the mass, $\gamma$ is the friction coefficient, $F$ is a constant force, and $\eta_{i}$’s are Gaussian noise satisfying $\left\langle\eta_{i}(0)\eta_{j}(t)\right\rangle=2m_{i}\gamma k_{B}T\delta(t)\delta_{ij}$. The probability of observing a given trajectory can be written in the Onsager-Machlup form Onsager and Machlup (1953)

[TABLE]

The relative stochastic action can be derived from the ratio of path probability with and without $F$,

[TABLE]

As only the velocities are asymmetric upon time reversal, we identify the time asymmetric and symmetric parts as

[TABLE]

As seen, the time symmetric part has a part proportional to $F$, and a part proportional to $F^{2}$; in other words, $S[\tilde{x}]=\sum_{i=1}^{p}t_{N}Q_{i}X^{i}$ where $p=2$ . Defining the affinity $X=F/2$, we arrive at the expressions of the dynamical variables,

[TABLE]

As before, $q_{2}(\tilde{x}_{t})=-\sum_{i=1}^{N}1/{m_{i}\gamma}$ is a constant independent of the dynamics. Therefore, the expression for nonequilibrium current can be written in a closed form

[TABLE]

where $\beta=1/k_{B}T$, and the second sum is taken over all combinations of nonnegative indices $i,j$ such that $i+j=n$. $\kappa_{i,j}^{n}$ denotes the $n$-th order cumulant, with the first few terms:

[TABLE]

Acknowledgements

The authors thank the UC Berkeley College of Chemistry for support. CYG was supported by the Kavli Energy NanoScience Institute.