Fixing the flux: a Dual Approach to Computing Transport Coefficients

Transport coefficients: a computational challenge

Formal setting: linear responses of stochastic dynamics

• Perturbed stochastic differential equation: \[ \mathrm{d}X^\eta_t = b(X^\eta_t)\,\mathrm{d}t+\sigma(X^\eta_t)\,\mathrm{d}W_t +\eta F(X^\eta_t)\,\mathrm{d}t,\qquad X^\eta_t\in\mathcal X. \] Parameters: vector fields \(\mathcal X\to T\mathcal X\), \(b\) (equilibrium drift), \(F\) (nonequilibrium forcing), \(\sigma\) (diffusion matrix), \(W\) Wiener process.

• Strength of the perturbation: \(\eta\). Equilibrium when \(\eta = 0\).

• Nonequilibrium ensemble: \(\psi_\eta\in\mathcal M_1(\mathcal X)=\) invariant measure for \(X^\eta\). Generally non-explicit for \(\eta\neq 0\).

• Define a response observable of flux \(R:\mathcal X \to\mathbb R\), with \(\mathbb E_{\psi_0} R = 0\)

• Linear response is the derivative of the expected flux at equilibrium: \[ \begin{equation}\alpha(F,R) = \underset{\eta \to 0}{\lim}\, \frac1\eta\mathbb E_{\psi_\eta}[R] \end{equation} \]

Classical estimators of the linear response

• Nonequilibrium static averages (NEMD): \[ \begin{equation}\widehat\alpha_{T,\eta}^{\mathrm{NEMD}}(F,R) = \frac1{T\eta}\int_0^T R(X_t^\eta)\,\mathrm{d}t,\quad \eta>0. \end{equation} \]
  
• Variance scales like \(T^{-1}\eta^{-2}\): simulation times \(\mathcal{O}(\eta^{-2}\varepsilon^{-2})\) for confidence intervals of width \(\mathcal{O}(\varepsilon)\).
• Equilibrium dynamical averages (Green-Kubo): \[ \begin{equation} \widehat\alpha_T^{\mathrm{GK}}(F,R) = \int_0^{T} \mathbb E_{\psi_0}[R(X_t^0)S(X_0^0)]\,\mathrm{d}t, \end{equation} \] where \(S\) is the so-called conjugate response, which is explicit in terms of \(F\).

Dual approach: fixing the flux

• Macroscopically: does the forcing cause a flux in the nonequilibrium steady-state, or does the flux cause a resisting forcing ?

• Microscopically: dynamics with fixed forcing / variable flux vs dynamics with a fixed flux / variable forcing.

• Constant-flux dynamics \[ \mathrm{d}Y^r_t = b(Y^r_t)\,\mathrm{d}t+\sigma(Y^r_t)\,\mathrm{d}W_t + F(Y^r_t)\mathrm{d}\Lambda^r_t,\quad \mathrm{d}\Lambda^r_t = \lambda(Y_t^r)\,\mathrm{d}t + \widetilde{\lambda}(Y_t^r)\,\mathrm{d}W_t. \]

• The forcing process/stochastic Lagrange multiplier \(\Lambda^r\) is constructed to enforce the constant flux condition: \(R(Y_t^{r}) = r \in {\mathbb{{R}}}\) for all \(t>0\).

• \(\Lambda^r\) satisfies a SDE with explicit coefficients \(\lambda,\widetilde{\lambda}\), which can be obtained under the controllability condition \(F^\top\nabla R\neq 0\).

• Generalization of the “Norton” method of Evans & Morriss (1980s) to stochastic dynamics.

• Also works for vector-valued fluxes/time-dependent constraints.

Dual approach: fixing the flux

Linear response

• Denote by \(\psi^r\) the constant-flux steady-state (or “Norton ensemble”).
• Dual definitions of the linear response \[ \alpha^*(F,R) = \underset{r\to 0}{\lim}\, \frac{r}{\mathbb E_{\psi^r}[\lambda]},\qquad \alpha(F,R) = \underset{\eta\to 0}{\lim}\,\frac{\mathbb E_{\psi_\eta}[R]}{\eta}. \]
• Set \(r(\eta):=\mathbb{E}_{\psi_\eta}[R],\,\eta(r):=\mathbb E_{\psi^r}[\lambda]\). If \(\eta(r(\eta)) = \eta + o(\eta)\) and \(r(\eta(r))=r + o(r)\), then \(\alpha^*(F,R)=\alpha(F,R)\).
• Reciprocal estimator of the linear response: \[ \begin{equation}\widehat\alpha_{T,r}^{*}(F,R) = \left(\frac1{T r}\int_0^T \lambda(Y_t^r)\,\mathrm{d}t\right)^{-1},\quad r>0. \end{equation} \]
• Explicit expression for \(\lambda,\widetilde{\lambda}\): \[ \lambda = -\frac1{F^\top \nabla R}\left(b^\top \nabla R +\frac12 \mathrm{Tr}\left[\sigma^\top \overline{P}_{F,\nabla R}^{\top}\nabla^2 R \overline{P}_{ F,\nabla R}\sigma\right]\right),\qquad \overline{P}_{F,\nabla R} = \mathrm{Id} - \frac{F\nabla R^\top}{F^\top\nabla R} \]
• \(\rightarrow\) link with variance reduction: uses \(-\frac1T\int_0^T \widetilde{\lambda}(Y_t^r)\mathrm{d}W_t\) as a control variate.

Example: mobility of a Langevin particle

• Specialize to underdamped Langevin dynamics. Coefficients \(b(q,p)=\begin{pmatrix} p \\-\nabla V(q) - \gamma p\end{pmatrix}\), \(\sigma(q,p) = \sqrt{2\gamma}\begin{pmatrix} 0 & 0 \\ 0 & \mathrm{Id} \end{pmatrix}\) in phase space \((q,p)\in ({\mathbb{{R}}}/\mathbb Z)^d\times {\mathbb{{R}}}^d\).
• Forcing/Flux pair: fixed direction \(\mathcal F = \begin{pmatrix}0\\ F\end{pmatrix} \in {\mathbb{{R}}}^{2d}\). Response is the system’s velocity in the \(F\)-direction: \(R(q,p) = F^\top p\). Controllability condition: \(\mathcal F^\top \nabla R = |F|^2>0\).
• Apply Itô’s formula to the constant-flux condition \(R(q^r_t,p^r_t) = r\): \[ \nabla R(q^r_t,p^r_t)^\top\left(\left[b(q_t^r,p_t^r) + \lambda(q_t^r,p_t^r)\mathcal F\right]\,\mathrm{d}t + \left[\sigma+\mathcal F\widetilde\lambda(q_t^r,p_t^r)\right]\mathrm{d}W_t\right) + \frac12\sigma\sigma^\top \nabla^2 R(q_t^r,p_t^r)\mathrm{d}t = 0 \]
• \[ \implies F^\top\left[-\nabla V(q^r_t) -\gamma p_t^r + \lambda(q_t^r,p_t^r)F\right] = 0,\qquad F^\top\left[\sqrt{2\gamma}\begin{pmatrix} 0 & 0 \\ 0 & \mathrm{Id} \end{pmatrix}+\mathcal F\widetilde\lambda(q_t^r,p_t^r)\right] = 0 \]
• \[ \implies \lambda(q,p) = -\frac{F^\top (-\nabla V(q)-\gamma p)}{|F|^2},\qquad \widetilde\lambda(q,p) = -\frac{\sqrt{2\gamma}}{|F|^2}\begin{pmatrix}0_q^\top & F^\top \end{pmatrix}. \]

Norton-Langevin dynamics for mobility

• Plugging these in, we get a dynamics \[ \left\{\begin{aligned} \mathrm{d}q_t^r &= p_t^r\,\mathrm{d}t\\ \mathrm{d}p_t^r &= \overline{P}_{F}\left(-\nabla V(q^r_t)\mathrm{d}t - \gamma p_t^r \mathrm{d}t + \sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_t\right) \end{aligned}\right.,\qquad \overline{P}_F = \mathrm{Id}_p - P_F,\qquad P_F = \frac{F F^\top}{F^\top F} \]
• Projection of the equilibrium dynamics on the constant-flux manifold.

General case: oblique projection

General case: oblique projection

General case: oblique projection

General case: oblique projection

General case: oblique projection

General case: oblique projection

Does it work ?

Does \(\alpha(F,R) = \alpha^*(F,R)\) ?

Does it always work ?

• No.
• Take \(F = \begin{pmatrix} 1 & 0 & 0 & 0 & \cdots\end{pmatrix}^\top\): “single-drift” forcing.
• \(\rightarrow\) question of sufficient conditions for the equivalence to hold.

Is it efficient ?

• Assuming equivalence of the constant-forcing and constant-flux ensembles, we have a choice to measure a given response property.
• Which ensemble to choose: the one with lower asymptotic variance for linear response estimators.

Some perspectives

• Many theoretical questions are open:
  • Ergodicity of the Norton dynamics
  • Equivalence of linear responses
  • Equivalence of nonequilibrium ensembles
  • Forcing fluctuations in the Norton ensemble, anomalous scaling
• Some ongoing work by Shiva Darshan & Gabriel Stoltz.
• Application to DPD: Xiaocheng Shang & Xinyi Wu (next talk!).