# Fixing the flux: A dual approach to computing transport coefficients

N. Blassel1,2 and G. Stoltz1,2
1: CERMICS, École des Ponts, Marne-la-Vallée, France
2: MATHERIALS project-team, Inria Paris, France

###### Abstract

We present a method to compute transport coefficients in molecular dynamics.
Transport coefficients quantify the linear dependencies of fluxes in non-
equilibrium systems subject to small external forcings. Whereas standard non-
equilibrium approaches fix the forcing and measure the average flux induced in
the system driven out of equilibrium, a dual philosophy consists in fixing the
value of the flux, and measuring the average magnitude of the forcing needed
to induce it. A deterministic version of this approach, named Norton dynamics,
was studied in the 1980s by Evans and Morris. In this work, we introduce a
stochastic version of this method, first developing a general formal theory
for a broad class of diffusion processes, and then specializing it to
underdamped Langevin dynamics, which are commonly used for molecular dynamics
simulations. We provide numerical evidence that the stochastic Norton method
provides an equivalent measure of the linear response, and in fact demonstrate
that this equivalence extends well beyond the linear response regime. This
work raises many intriguing questions, both from the theoretical and the
numerical perspectives.

##  1 Introduction

Molecular dynamics (MD), much like computational statistical physics at large,
aims at predicting macroscopic properties of a molecular medium by the use of
computer simulations. In order to do so, a classical model of the interactions
between atoms is constructed, and dynamical evolution laws are specified. One
can then simulate typical trajectories of the system through phase space,
which, provided they are long enough, allows to estimate to a prescribed level
of statistical accuracy the quantities of interest. See [3] for a general
theoretical account of statistical physics, [45, 1, 16] for an emphasis on
numerical methods, and [29] for an overview of the mathematical aspects of
molecular simulation. MD simulations may prove useful when the properties of
interest are impractical to measure due to physical or cost constraints
associated with the experimental setup, or alternatively serve as surrogate
tests for theoretical models. At any rate, MD has grown to occupy an important
role in many applications ranging from pharmacology and molecular biology to
materials science and condensed matter physics, besides having had a
significant impact on statistical methodology at large, through the
dissemination of tools such as the Metropolis–Hastings algorithm [35]. First
applications of MD included the computation of static properties at
equilibrium, in particular thermodynamic quantities or free energy
differences, which still are of considerable interest today. We refer to [7]
for a historical perspective.

Another, more difficult problem is the measurement of dynamical properties,
quantities which depend on the trajectory itself, and capture the behavior of
systems undergoing a macroscopic time evolution: this is the object of non-
equilibrium statistical physics. Of particular interest is the computation of
transport coefficients, giving a measure of the sensitivity in the response of
a system at equilibrium to the magnitude of perturbations driving it out of
equilibrium. See [14, 44] for overviews, and [34, Section 5] for a
mathematical presentation. Transport coefficients are of practical importance,
since they characterize the diffusion, heat conduction or viscosity properties
of a molecular medium, and enter as parameters in macroscopic evolution
equations such as the Navier–Stokes equation. One standard approach to compute
transport coefficients by molecular dynamics relies on the celebrated
Green–Kubo formula [18], which expresses transport coefficients as integrated
time-correlations between appropriate fluxes in the system at equilibrium.
Another standard approach, see for instance [6], is to directly simulate the
non-equilibrium perturbation, and to measure the resulting average response,
which is, at the macroscopic level and in the limit of a small perturbation,
proportional to the magnitude of this perturbation: this is the so-called non-
equilibrium molecular dynamics (NEMD) approach, reviewed in [7]. Estimators
deriving from these approaches however suffer from large statistical errors,
as quantified in [41, Proposition 2.4] for instance, and convergence requires
the simulation of very long trajectories, which comes at a high computational
cost. A key metric to measure this cost is, for any given method, the
asymptotic variance of estimators of the transport coefficient. Although some
variance reduction techniques have been proposed to compute transport
coefficients (see [42] for a recent review), efficiently estimating these
quantities is are still an important area of research.

From a thermodynamic point of view, the NEMD approach can be understood as
fixing the magnitude of the non-equilibrium forcing, and measuring the
resulting flux in the system driven out of equilibrium, which is conceived as
a microscopic state variable. For a small enough magnitude of the forcing,
this flux responds approximately linearly to the forcing, and the transport
coefficient is precisely the associated proportionality constant. This method
also yields a construction of the non-equilibrium ensemble, by defining it as
the steady state of the system evolving according to the microscopic dynamics.

However convenient from an implementation standpoint, there is no physical
reason to a priori favor an ensemble in which the forcing field is exactly
fixed. At the macroscopic level, fluxes and forces play a symmetric role,
which opens the possibility for dual approaches in the computation of non-
equilibrium responses. Other approaches have been proposed to construct non-
equilibrium ensembles in which the forcing does not play such a distinguished
role. One can cite McLennan ensembles [5], expressing non-equilibrium steady-
states as first-order perturbations of the Boltzmann distribution, where the
leading order correction term is the time-integral of some conjugate response;
or the dynamical approaches by Komorowski and Olla [26], whereby an average
flux is induced by a time-periodic forcing. When considering a dual
perspective, the most radical approach is to fix the flux exactly, and measure
the average magnitude of the forcing needed to induce it in the microscopic
dynamics. This yields an alternative measure of the transport coefficient. By
analogy with Ohm’s law, one can think of the NEMD approach as a way to measure
a conductance, and the radical dual approach as a way to measure a resistance.

The microscopic implementation of this dual approach is what we will refer to
as the Norton method, in accordance with the terminology proposed by Evans and
Morriss in [15], by reference to the Norton–Thévenin reciprocity from
electrical circuit theory. This approach amounts to considering a constrained
dynamics on a submanifold of phase space consisting of a level-set of the flux
observable. The idea of using constrained dynamics to simulate non-equilibrium
systems already appeared in the 1983 work [13], where it was applied to
capture the mobility, and was also explored in [22]. It was further applied to
shear viscosity computations in [12], where the consistency of the approach
was demonstrated. From a theoretical perspective, formal results under ergodic
hypotheses were obtained, including linear response expressions for the
transport coefficient in [15], as well as a result on the equivalence of non-
equilibrium ensembles in [11].

However, Norton dynamics were only considered in a deterministic setting, and,
despite promising results, their potential for practical use was not fully
explored. In particular, the numerical performance of estimators of transport
coefficients based on time averages of this dynamics has not, to the best of
our knowledge, been studied yet. Our aim is to extend the Norton philosophy to
the stochastic setting, both for academic motivations (obtaining new results
on equivalence of ensembles in non-equilibrium systems), and for numerical
reasons, as Norton dynamics potentially allow to more efficiently compute
transport coefficients. We will consider both general diffusion processes, and
Langevin-type dynamics, which are commonly used in molecular dynamics [36].

#### Contributions of this work.

Let us highlight our main contributions:

  * $\circ$

We define a stochastic version of the Norton method, which is convenient from
the theoretical point of view as rigorous ergodicity results exist for
stochastic dynamics, and are also relevant for simulation since stochastic
dynamics are nowadays very commonly used in molecular dynamics.

  * $\circ$

We specialize the Norton method to underdamped Langevin dynamics, and apply it
to compute the mobility and shear viscosity of a fluid. We demonstrate
numerically that the Norton method gives consistent estimations of the linear
response in the thermodynamic limit.

  * $\circ$

We observe on the numerical examples we consider that the non-equilibrium
responses coincide in fact far outside of the linear regime, raising the
question of equivalence between non-equilibrium ensembles.

  * $\circ$

We offer numerical evidence that, in some situations, the Norton method is
preferable to the standard NEMD approach, in the sense that the Norton
estimators for the transport coefficients lead to estimates with a smaller
statistical error than their NEMD counterparts.

Many points in the mathematical analysis of the Norton approach are left open
at this stage. This should be seen as an invitation to further study the
properties of these intriguing dynamics

#### Outline.

This paper is organized as follows. We recall in Section 2 a general framework
for stochastic dynamics out of equilibrium. In Section 3, we introduce the
Norton approach, deriving an expression for the diffusion process defining the
dynamics, before generalizing the approach to the case of multiple constraints
and time-dependent fluxes. In Section 4, we specialize the setting to non-
equilibrium kinetic Langevin dynamics used for mobility and shear viscosity
computations, describing how to apply the Norton philosophy. We also give a
physical interpretation of the Norton dynamics, as one satisfying an oblique
version of Gauss’s principle of least constraint. In Section 5, we discuss how
to discretize stochastic Norton dynamics, describing in particular a family of
schemes obtained by an operator splitting approach. We present the results of
our numerical experiments in Section 6, demonstrating the consistency with
usual NEMD in the thermodynamic limit, and motivating that Norton simulations
can be more efficient than NEMD ones. We conclude in Section 7 by discussing
the many open questions and perspectives raised by the Norton method in the
stochastic context.

##  2 Non-equilibrium molecular dynamics

We recall in this section the standard framework of NEMD, first presenting the
reference dynamics, before introducing the non-equilibrium dynamics and
defining transport coefficients. We finally discuss the statistical properties
of NEMD estimators of the transport coefficient. We refer the interested
reader to [34, Section 5] and [41] for a more in-depth discussion of the
mathematical properties of the NEMD method.

#### Reference dynamics.

We consider a class of non-equilibrium ensembles, obtained as the steady-
states of certain stochastic processes corresponding to a perturbation of a
reference process (usually a dynamics at equilibrium). We denote by
$\mathcal{X}$ the state space of the system. Typically, we consider
$\mathcal{X}=\mathbb{R}^{d}$ or $\mathbb{T}^{d}$, with
$\mathbb{T}=\mathbb{R}/\mathbb{Z}$ the one-dimensional torus, or a product of
the two. We introduce a smooth vector field $b:\mathcal{X}\to\mathbb{R}^{d}$,
corresponding to the equilibrium drift, and a matrix valued map
$\sigma:\mathcal{X}\to\mathbb{R}^{d\times m}$ corresponding to the diffusion.
The reference dynamics is the following stochastic differential equation
(SDE):

$\mathrm{d}X_{t}=b(X_{t})\,\mathrm{d}t+\sigma(X_{t})\,\mathrm{d}W_{t},$ (1)

where $(W_{t})_{t\geqslant 0}$ a standard $m$-dimensional Brownian motion.
Common choices include the overdamped Langevin dynamics

$\mathrm{d}X_{t}=-\nabla
V(X_{t})\,\mathrm{d}t+\sqrt{\frac{2}{\beta}}\,\mathrm{d}W_{t},$

and the underdamped or kinetic Langevin dynamics (see Equation (35) below).

#### Non-equilibrium perturbations.

We consider in this work the case where the perturbation arises from an
external non-gradient forcing in the drift of the underlying diffusion
process, determined by a smooth vector field $F:\mathcal{X}\to\mathbb{R}^{d}$.
The non-equilibrium dynamics is given by the following SDE:

$\mathrm{d}X^{\eta}_{t}=(b+\eta
F)(X^{\eta}_{t})\,\mathrm{d}t+\sigma(X^{\eta}_{t})\,\mathrm{d}W_{t}.$ (2)

The parameter $\eta\in\mathbb{R}$ is a scalar modulating the strength of the
perturbation. The equilibrium dynamics is recovered in the absence of a non-
equilibrium forcing, i.e. $\eta=0$. The infinitesimal generator of the
dynamics (2) can be decomposed as the sum

$\mathcal{L}_{\eta}=\mathcal{L}+\eta\widetilde{\mathcal{L}},\qquad\mathcal{L}=b\cdot\nabla+\frac{1}{2}\sigma\sigma^{\intercal}:\nabla^{2},\qquad\widetilde{\mathcal{L}}=F\cdot\nabla,$
(3)

where $\nabla^{2}$ denotes the Hessian matrix and $:$ denotes the Frobenius
inner product $A:B=\mathrm{Tr}(A^{\intercal}B)$. Note that $\mathcal{L}$ is
the generator of the reference dynamics, and $\widetilde{\mathcal{L}}$ encodes
its perturbation. The invariant probability measure satisfies the stationary
Fokker-Planck equation

$\mathcal{L}_{\eta}^{\dagger}\psi_{\eta}=0,$ (4)

where $\mathcal{L}_{\eta}^{\dagger}$ is the flat $L^{2}(\mathcal{X})$-adjoint
of the generator. Existence and regularity results for solutions of (4) depend
on the particular form of the dynamics, as do properties pertaining to
convergence to the steady-state. It is often possible to leverage the standard
analytical framework of parabolic and elliptic partial differential equations
to prove this convergence, although for degenerate diffusions such as the
underdamped Langevin dynamics, one has to resort to more sophisticated tools,
such as hypoellipticity [23] to obtain regularity of the steady-state, or
hypocoercivity [47] to show exponential decay of the evolution semigroup.

We assume in the remainder of this section that, for $|\eta|$ small enough,
the dynamics admits a unique invariant probability measure for which it is
ergodic, and denote the corresponding expectation by $\mathbb{E}_{\eta}$.
Given a response observable $R:\mathcal{X}\to\mathbb{R}$ such that

$\mathbb{E}_{0}\left[R\right]=0,$ (5)

which we think of as measuring a flux in the system out of equilibrium, we
define the associated transport coefficient as the following limit, provided
it is well defined:

$\alpha=\underset{\eta\to
0}{\lim}\frac{\mathbb{E}_{\eta}\left[R\right]}{\eta}.$ (6)

Rigorous assumptions under which this limit exists are given in [21]. This
definition suggests a simple and natural method to estimate these
coefficients: one can compute ergodic averages of $R$ over trajectories of the
non-equilibrium dynamics (2), and estimate the linear relation between $\eta$
and $R$ for one or several values of $\eta$ in the linear response regime. The
finite-difference estimator for the limit (6) is given by the following
ergodic average:

$\widehat{\alpha}_{T,\eta}=\frac{1}{\eta
T}\int_{0}^{T}R(X^{\eta}_{t})\,\mathrm{d}t.$ (7)

The consistency of such estimators is a consequence of the pathwise ergodicity
for the process (2), which in many cases can be proven using the results of
[25]. The latter result also implies the uniqueness of the steady-state. The
existence of the steady-state is often obtained by Lyapunov techniques, see
for instance [37, Theorem 8.3].

#### Statistical properties of the estimator $\widehat{\alpha}_{T,\eta}$.

A challenge posed by the NEMD method is that the estimator (7) is plagued by
large statistical errors when $|\eta|$ is small, which is often required to
remain in the linear response regime, and avoid biases arising from nonlinear
terms in the response. More precisely, one can show that the asymptotic
variance of the estimator (7) scales as $\eta^{-2}$ as $\eta$ approaches $0$.
Indeed, under technical decay conditions on the evolution semigroup generated
by (3), one can easily show that the asymptotic variance is given by

$\sigma_{\eta}^{2}=\lim_{T\to\infty}\,T\mathrm{Var}_{\eta}\left(\widehat{\alpha}_{T,\eta}\right)=\frac{2}{\eta^{2}}\int_{0}^{\infty}\mathbb{E}_{\eta}[R(X_{t}^{\eta})R(X_{0}^{\eta})]\,\mathrm{d}t,$
(8)

where $\mathrm{Var}_{\eta}$ denotes the variance with respect to
$\mathbb{E}_{\eta}$. Defining the correlation time by

$\Theta_{\eta}(R)=\int_{0}^{\infty}\frac{\mathbb{E}_{\eta}[R(X_{t}^{\eta})R(X_{0}^{\eta})]}{\mathbb{E}_{\eta}[R^{2}]}\,\mathrm{d}t,$
(9)

we further get, using a first-order expansion in powers of $\eta$ (whose
validity has to be verified on a case-by-case basis),

$\sigma_{\alpha,\eta}^{2}=\frac{2}{\eta^{2}}\mathrm{Var}_{\eta}(R)\Theta_{\eta}(R)=\frac{2}{\eta^{2}}\mathrm{Var}_{0}(R)\Theta_{0}(R)+\mathrm{O}\left(\frac{1}{\eta}\right).$
(10)

In summary, the asymptotic variance is, at dominant order in $|\eta|$, the
asymptotic variance of the time averages of $R$ under the equlibrium dynamics,
divided by $\eta^{2}$.

The leading order of the asymptotic variance highlights why it is often
computationally expensive to obtain accurate estimates of transport
coefficients. Although the computation of such coefficients is recognized as
being a difficult problem in practice, only a handful of research works have
been proposed to alleviate this issue, including replacing the forcing $F$ by
a so-called synthetic forcing devised to induce the same transport coefficient
while the range of the linear response [41], or relying on carefully
constructed couplings [9]. We refer the reader to [42, Section 4] for a recent
overview of current variance reduction techniques in non-equilibrium molecular
dynamics.

##  3 A stochastic Norton method

The Norton approach exploits the macroscopic duality between thermodynamic
forces and fluxes: at the macroscopic level, one can equivalently choose to
measure the current induced by a constant force, or the resistance opposed to
a constant current. The microscopic translation of this duality is the
introduction of a new non-equilibrium ensemble in which the flux is held
fixed. As in the NEMD case, we define this ensemble as the invariant
probability measure for a particular stochastic process, which we refer to as
the Norton dynamics.

In Section 3.1, we present the Norton pertubation approach for a generic
reference dynamics of the form (2). We then proceed in Section 3.2 to write
the dynamics in closed form, by making explicit the constraining force on the
flux. In Section 3.3, we give the expression of the Norton analogs of the
transport coefficient (6), and discuss how their statistical properties can be
formally analyzed. We finally show in Section 3.4 how the Norton approach can
be extended to the cases of multiple forcings, or time-dependent flux
constraints.

###  3.1 Presentation of the dynamics

At the dynamical level, the Norton ensemble is defined as the invariant
probability measure of the following stochastic differential equation:

$\left\\{\begin{aligned}
\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},\\\
R(Y^{r}_{t})&=R(Y^{r}_{0})=r.\end{aligned}\right.$ (11)

Here, the evolution of the state is given by the dynamics of
$Y^{r}_{t}\in\mathcal{X}$, and $r\in\mathbb{R}$ is the magnitude of the
response flux, which is maintained constant. The stochastic dynamics therefore
evolves on the submanifold

$\Sigma_{r}=\left\\{y\in\mathcal{X},\ R(y)=r\right\\}=R^{-1}\\{r\\}$ (12)

of the full state space. The dynamics (11) can formally still be considered as
a perturbation of the equilibrium dynamics, in the same direction as the
Thévenin process (2), but with a random intensity given by the stochastic
process $\Lambda^{r}_{t}$, acting as the control by which the constant-flux
condition is enforced. Provided $\mathbb{E}_{0}[R]=0$, we can further
interpret the Norton dynamics $Y_{t}^{0}\in\Sigma_{0}$ as an equilibrium
dynamics, constrained to exactly preserve the flux. The relationship between
the equilibrium ensemble in which the response fluctuates, and the Norton
equilibrium ensemble in which it is exactly fixed at zero, is reminiscent of
the relationship between the canonical and microcanonical ensembles when $R$
is a spatial average over local quantities.

We next to show that $\Lambda^{r}$ can in fact be defined as an Itô process
adapted to the natural filtration of the $m$-dimensional Brownian motion $W$.
More precisely, we decompose the intensity of the forcing using the following
ansatz:

$\Lambda^{r}_{t}=\Lambda^{r}_{0}+\int_{0}^{t}\lambda(Y^{r}_{s})\,\mathrm{d}s+\widetilde{\Lambda}_{t}^{r},\qquad\widetilde{\Lambda}_{t}^{r}=\int_{0}^{t}\widetilde{\lambda}(Y^{r}_{s})\,\mathrm{d}W_{s},$
(13)

where $\lambda,\widetilde{\lambda}$ are functions defined on $\mathcal{X}$,
with $\lambda$ taking values in $\mathbb{R}$ and $\widetilde{\lambda}$ in
$\mathbb{R}^{1\times m}$. This ansatz is natural for constrained dynamics (see
for instance [32, Chapter 3]), and will be confirmed a posteriori (see in
particular (15) and (20) in Section 3.2 below). The average forcing in the
Norton ensemble can then be defined as the expectation of $\lambda$ under the
steady-state probability measure, neglecting the zero-mean contribution of
$\widetilde{\lambda}$ (see Section 3.3). Numerically, these averages can be
computed as ergodic averages over discretized trajectories of the Norton
dynamics, as discussed in Section 5.3.

###  3.2 A closed form for the forcing process

We make precise here the expressions of the function $\lambda(Y^{r}_{t})$ and
the martingale $\widetilde{\Lambda}^{r}_{t}$ in (13), which allows us to write
the Norton dynamics without explicit reference to the forcing. We assume that
$\Lambda^{r}$ is of the form (13), and verify a posteriori that this ansatz is
valid. Applying Itô’s formula to the constant response condition
$R(Y^{r}_{t})=r$ yields

$\nabla
R(Y^{r}_{t})\cdot\left[b(Y^{r}_{t})\,\mathrm{d}t+\sigma(Y^{r}_{t})\,\mathrm{d}W_{t}+F(Y^{r}_{t})\,\mathrm{d}\Lambda^{r}_{t}\right]+\frac{1}{2}\nabla^{2}R(Y^{r}_{t}):\mathrm{d}\left\langle
M^{r}\right\rangle_{t}=0,$ (14)

where $\langle M\rangle_{t}$ denotes the quadratic covariation process for the
martingale part in the Itô decomposition of $Y^{r}$:

$\mathrm{d}M^{r}_{t}=\sigma(Y^{r}_{t})\,\mathrm{d}W_{t}+F(Y^{r}_{t})\,\widetilde{\lambda}(Y^{r}_{t})\,\mathrm{d}W_{t}.$
(15)

Using the uniqueness of the Itô decomposition, we can identify the martingale
increment in (14) as

$\mathrm{d}\widetilde{\Lambda}^{r}_{t}=-\frac{\nabla
R(Y^{r}_{t})\cdot\sigma(Y^{r}_{t})\,\mathrm{d}W_{t}}{\nabla R(Y^{r}_{t})\cdot
F(Y^{r}_{t})},$ (16)

provided that $\nabla R\cdot F\neq 0$ almost surely, which we assume here and
in the sequel. Plugging this equality in (15) in turn gives

$M_{t}^{r}=\int_{0}^{t}\left(\mathrm{Id}-\frac{F(Y^{r}_{s})\otimes\nabla
R(Y^{r}_{s})}{F(Y^{r}_{s})\cdot\nabla
R(Y^{r}_{s})}\right)\sigma(Y^{r}_{s})\,\mathrm{d}W_{s}=\int_{0}^{t}\overline{P}_{F,\nabla
R}\sigma(Y^{r}_{s})\,\mathrm{d}W_{s},$

so that the covariation of the martingale increment is

$d\left\langle M^{r}\right\rangle_{t}=\left[\overline{P}_{F,\nabla
R}\sigma\sigma^{\intercal}\overline{P}_{\nabla
F,R}^{\intercal}\right](Y^{r}_{t})\,\mathrm{d}t.$ (17)

In the latter two expressions, we make use of the following non-orthogonal
projector-valued maps, defined for vector fields $A,B$ such that $A(x)\cdot
B(x)\neq 0$:

$P_{A,B}(x)=\frac{A(x)\otimes B(x)}{A(x)\cdot
B(x)},\qquad\overline{P}_{A,B}(x)=\mathrm{Id}-P_{A,B}(x).$ (18)

The action of the projector is given, for $\xi\in\mathbb{R}^{d}$ and
$x\in\mathcal{X}$, by

$P_{A,B}(x)(\xi)=\frac{B(x)\cdot\xi}{A(x)\cdot B(x)}A(x).$

For notational convenience, we introduce

$\Pi_{F,\nabla R,\sigma}(y)=\left[\overline{P}_{F,\nabla
R}\sigma\sigma^{\intercal}\overline{P}_{F,\nabla R}^{\intercal}\right](y).$
(19)

We next proceed to identifying the bounded-variation increments on both sides
of (14). After rearrangement and substitution of (17), one obtains the
following expression:

$\lambda_{t}^{r}=\lambda(Y_{t}^{r}),\qquad\lambda=-\frac{1}{F\cdot\nabla
R}\left(b\cdot\nabla R+\frac{1}{2}\nabla^{2}R:\Pi_{F,\nabla R,\sigma}\right).$
(20)

Substituting the expression for $\Lambda^{r}_{t}$ in (11) yields the following
expression for the Norton dynamics:

$\mathrm{d}Y^{r}_{t}=\overline{P}_{F,\nabla
R}(Y^{r}_{t})\left[b(Y^{r}_{t})\mathrm{d}t+\sigma(Y^{r}_{t})\,\mathrm{d}W_{t}\right]-\frac{\left(\nabla^{2}R:\Pi_{F,\nabla
R,\sigma}\right)(Y^{r}_{t})}{2\nabla R(Y^{r}_{t})\cdot
F(Y^{r}_{t})}F(Y^{r}_{t})\,\mathrm{d}t.$ (21)

It can now be checked a posteriori, using Itô’s formula, that the dynamics
(21) is such that $R(Y^{r}_{t})=r$ for all $t\geqslant 0$, provided the
coefficients are smooth. In Figure 1, we illustrate geometrically the action
of the projector $\overline{P}_{F,\nabla R}$ at a point $y$ on a vector
$\Delta y$.

$y$$\Sigma_{r}$$T_{\Sigma_{r}}(y)$$\nabla R(y)$$F(y)$$\Delta
y$$\overline{P}_{F,\nabla R}(y)\Delta y$ Figure 1: Action of the projector
$\overline{P}_{F,\nabla R}$: the increment $\Delta y$ is projected onto the
tangent space $T_{\Sigma_{r}}(y)=\\{z\in\mathbb{R}^{d}\,|\,\nabla R(y)\cdot
z=0\\}$ in the direction $F(y)$.

Without further specifying the particular choice for the reference dynamics,
the response flux observable $R$ and non-equilibrium forcing $F$, it is
difficult to make general comments about the well-posedness of (21). Let us
however emphasize that a crucial condition for the dynamics to be well-defined
is that the denominator $F(Y^{r}_{t})\cdot\nabla R(Y^{r}_{t})$ in the
expression of the projector $P_{F,\nabla R}$ and in the last term of (21)
(which represents some curvature correction) should not vanish. Thinking of
the extreme case where $\nabla R$ and $F$ are everywhere orthogonal, we see
that this requirement translates into a controllability condition: in this
case, any forcing in the direction $F$ has no effect on the flux, thus there
is no way to control the latter using such a perturbation. More generally,
starting from a configuration for which $F\cdot\nabla R(q)=0$, it is not
possible to maintain the value of the response function using the forcing $F$.
We therefore assume in the sequel that the condition

$\forall y\in\Sigma_{r},\qquad F(y)\cdot\nabla R(y)\neq 0,$ (22)

is satisfied. We discuss the validity of this condition for the numerical
examples we consider in Section 6.

###  3.3 Norton analogs of the transport coefficient

We assume the well-posedness of the dynamics (11) (or equivalently (21)), and
also the existence and uniqueness of the invariant steady-state probability
measure for this dynamics, whose expectation is denoted by
$\mathbb{E}_{r}^{*}$. We assume that

$\mathbb{E}_{0}^{*}[\lambda]=0.$ (23)

This condition is the Norton counterpart of the centering condition (5) for
the observable in usual NEMD simulations. When (23) holds, the transport
coefficient for the Norton dynamics is defined by analogy with (6) as

$\alpha_{F,R}^{*}=\underset{r\to
0}{\lim}\,\frac{r}{\displaystyle{\mathbb{E}_{r}^{*}\left[\lambda\right]}},$
(24)

provided the limit exists. In equation (24), $\alpha_{F,R}^{*}$ can be
interpreted as the inverse of the resistance to the non-equilibrium forcing.
Note that the average forcing in the denominator of the right-hand side of
(24) only involves the bounded varation part of (13), since the expectation of
the martingale part vanishes.

Provided the Norton dynamics is ergodic with respect to the steady-state, a
natural estimator of the transport coefficient can be constructed by replacing
ensemble averages by trajectory averages, similarly to what is done in
standard NEMD simulations. More precisely, the estimator is computed using
ergodic averages of $\lambda(Y^{r}_{t})$ under $\mathbb{E}_{r}^{*}$, that is,

$\widehat{\alpha}_{T,r}^{*}=\frac{rT}{\displaystyle{\int_{0}^{T}\lambda(Y_{t}^{r})\,\mathrm{d}t}}.$
(25)

The statistical properties of this estimator can be analyzed similarly to what
is done for the estimator (7), leading to a result similar to (10). Assuming
that a central limit theorem holds, and that $F$, $R$ are such that
$\mathbb{E}_{0}^{*}[\lambda]=0$, we can formally write the asymptotic variance
associated with the estimator (25) as

$\displaystyle\sigma_{\alpha,r}^{2,*}$
$\displaystyle=\lim_{T\to\infty}T\mathrm{Var}^{*}_{r}(\widehat{\alpha}_{T,r}^{*})=\frac{r^{4}}{\mathbb{E}_{r}^{*}[\lambda]^{4}}\lim_{T\to\infty}T\mathrm{Var}^{*}_{r}\left(\frac{1}{\widehat{\alpha}_{T,r}^{*}}\right)$
(26)
$\displaystyle=\frac{r^{4}}{\mathbb{E}_{r}^{*}[\lambda]^{4}}\frac{2}{r^{2}}\mathrm{Var}_{r}^{*}(\lambda)\Theta_{r}^{*}(\lambda),$

where $\mathrm{Var}_{r}^{*}$ denotes the variance with respect to
$\mathbb{E}_{r}^{*}$, and the Norton correlation time is defined similarly to
the NEMD case by

$\Theta_{r}^{*}(\lambda)=\int_{0}^{\infty}\frac{\displaystyle{\mathbb{E}_{r}^{*}[\lambda(Y_{t}^{r})\lambda(Y_{0}^{r})]}}{\displaystyle{\mathbb{E}_{r}^{*}[\lambda^{2}]}}\,\mathrm{d}t.$
(27)

The second equality in (26) follows from the Delta method [46, Chapter 3]
applied to the reciprocal of $\widehat{\alpha}_{T,r}^{*}$. Using the first-
order expansion given by (24), and assumptions similar to the ones leading to
(10), we may further write, for $|r|$ small:

$\displaystyle\sigma_{\alpha,r}^{2,*}$
$\displaystyle=\frac{2\left(\alpha_{F,R}^{*}\right)^{4}}{r^{2}}\mathrm{Var}_{r}^{*}(\lambda)\Theta_{r}^{*}(\lambda)+\mathrm{O}\left(\frac{1}{r}\right)$
$\displaystyle=\frac{2\left(\alpha_{F,R}^{*}\right)^{4}}{r^{2}}\mathrm{Var}_{0}^{*}(\lambda)\Theta_{0}^{*}(\lambda)+\mathrm{O}\left(\frac{1}{r}\right),$

As in the NEMD setting, one has to verify the validity of each of the
expansion in powers of $|r|$ for each of the quantities considered. However,
whereas in the NEMD setting this may be done in some cases by using a
perturbative expansion of the non-equilibrium measure in powers of $\eta$ (as
done in [41], for instance), there is in the Norton setting an additional
technical difficulty in that the invariant measure for the Norton dynamics
(11) is supported on the $(d-1)$-dimensional manifold $\Sigma_{r}$, which is
disjoint from the Norton equilibrium manifold $\Sigma_{0}$, and also of zero
Lebesgue measure, so that it is both singular with respect to the invariant
measure of the Norton equilibrium dynamics and to that of the reference
dynamics (2). A dedicated analysis has therefore to be performed to justify
these approximations, and rigorously establish the validity of (non-)linear
response formulas.

###  3.4 Two straightforward generalizations

For ease of presentation, we restricted ourselves in Sections 3.1 to 3.3 to
the case where only one flux is fixed. However, the derivation of Norton
dynamics can be straightforwardly extended to the case where several fluxes
are simultaneously constrained. In fact, one can even consider situations
where the response depends on time. These two generalizations can of course be
combined.

#### Multiple constraints.

We consider in this paragraph the case where $c\geqslant 1$ forces
$F_{1},\dots,F_{c}$ act on the reference system with magnitudes
$\Lambda_{1,t}^{\mathbf{r}},\dots,\Lambda_{c,t}^{\mathbf{r}}$, chosen to
maintain constant the value of $c$ fluxes $R_{1},\dots,R_{c}$. The dynamics is
given by

$\left\\{\begin{aligned}
\mathrm{d}Y_{t}^{\mathbf{r}}&=b(Y_{t}^{\mathbf{r}})\,\mathrm{d}t+\sigma(Y_{t}^{\mathbf{r}})\,\mathrm{d}W_{t}+\sum_{i=1}^{c}F(Y_{t}^{\mathbf{r}})\,\mathrm{d}\Lambda_{i}^{\mathbf{r}},\\\
R_{i}(Y_{t}^{\mathbf{r}})&=r_{i},\qquad 1\leqslant i\leqslant
c,\end{aligned}\right.$

where $r_{i}$ gives the prescribed value of the corresponding flux $R_{i}$. By
defining the following maps defined on $\mathcal{X}$

$R=\begin{pmatrix}R_{1}\\\ \vdots\\\
R_{c}\end{pmatrix}\in\mathbb{R}^{c},\qquad F=\begin{pmatrix}|&&|\\\
F_{1}&\dots&F_{c}\\\ |&&|\end{pmatrix}\in\mathbb{R}^{d\times c},$

as well as the following vectors,

$\Lambda_{t}^{\mathbf{r}}=\begin{pmatrix}\Lambda_{1,t}^{\mathbf{r}}\\\
\vdots\\\
\Lambda_{c,t}^{\mathbf{r}}\end{pmatrix}\in\mathbb{R}^{c},\qquad\mathbf{r}=\begin{pmatrix}r_{1}\\\
\vdots\\\ r_{c}\end{pmatrix}\in\mathbb{R}^{c},$

the Norton dynamics can be written similarly to (11), upon replacing $r$ by
$\mathbf{r}$:

$\left\\{\begin{aligned}
\mathrm{d}Y^{\mathbf{r}}_{t}&=b(Y^{\mathbf{r}}_{t})\,\mathrm{d}t+\sigma(Y^{\mathbf{r}}_{t})\,\mathrm{d}W_{t}+F(Y^{\mathbf{r}}_{t})\,\mathrm{d}\Lambda^{\mathbf{r}}_{t},\\\
R(Y^{\mathbf{r}}_{t})&=R(Y^{\mathbf{r}}_{0})=\mathbf{r}.\end{aligned}\right.$
(28)

Computations identical to the ones leading to (21) may be performed. As these
are verbatim the same, we simply state the result and refer the reader to
Appendix A.1 for the complete derivation. The dynamics (28) can be written in
closed form, using non-orthogonal projectors, as

$\displaystyle\mathrm{d}Y^{\mathbf{r}}_{t}$
$\displaystyle=\overline{P}_{F,\nabla
R}(Y^{\mathbf{r}}_{t})\left[b(Y^{\mathbf{r}}_{t})\,\mathrm{d}t+\sigma(Y^{\mathbf{r}}_{t})\,\mathrm{d}W_{t}\right]$
(29) $\displaystyle\qquad-\frac{F(Y^{\mathbf{r}}_{t})}{2}(\nabla
R(Y^{\mathbf{r}}_{t})^{\intercal}F(Y^{\mathbf{r}}_{t}))^{-1}\left(\nabla^{2}R(Y^{\mathbf{r}}_{t})^{\intercal}:\Pi_{F,\nabla
R,\sigma}(Y^{\mathbf{r}}_{t})\right),$

where $\nabla R:\mathcal{X}\to\mathbb{R}^{d\times c}$ and
$\nabla^{2}R:\mathcal{X}\to\mathbb{R}^{d\times d\times c}$ are respectively
the Jacobian and Hessian matrices of the fluxes. Here, we define the
contraction product by

$\forall\,(A,B)\in\mathbb{R}^{d\times d\times c}\times\mathbb{R}^{d\times
d},\qquad A:B=\left(\sum_{j,k=1}^{d}A_{kji}B_{kj}\right)_{1\leqslant
i\leqslant c}\in\mathbb{R}^{c},$

the projector $\overline{P}_{F,\nabla R}$ is given, for
$A,B\in\mathbb{R}^{d\times c}$, by

$P_{A,B}=A(B^{\intercal}A)^{-1}B^{\intercal},\qquad\overline{P}_{A,B}=\operatorname{Id}-P_{A,B},$

and $\Pi_{F,\nabla R,\sigma}$ is defined similarly to (19) by

$\Pi_{F,\nabla R,\sigma}=\overline{P}_{F,\nabla
R}\sigma\sigma^{\intercal}\overline{P}_{F,\nabla R}^{\intercal}.$

Note that the invariant measure is supported on a codimension $c$ submanifold
$\Sigma_{\mathbf{r}}\subseteq\mathcal{X}$, and that the controllability
condition becomes

$\forall y\in\Sigma_{\mathbf{r}},\qquad\operatorname{det}(\nabla
R(y)^{\intercal}F(y))\neq 0.$ (30)

The average value of the forcing can be written as the ergodic average of the
following vector-valued observable, which is the analog of (20):

$\lambda=-(\nabla R^{\intercal}F)^{-1}\left[\nabla
R^{\intercal}b+\frac{1}{2}\left(\nabla^{2}R^{\intercal}:\Pi_{F,\nabla
R,\sigma}\right)\right].$ (31)

Dynamics such as (28) should in particular allow for the numerical computation
of Onsager relations in the Norton ensemble.

#### Time-dependent fluxes.

One can also extend the dynamics to the case where we replace the condition
$R(Y_{t}^{r})=r$ by $R(Y_{t}^{\mathcal{R}})=\mathcal{R}_{t}$, when
$\mathcal{R}_{t}$ is the Itô process defined by

$\mathcal{R}_{t}=R(Y^{\mathcal{R}}_{0})+\int_{0}^{t}\bar{r}_{s}\,\mathrm{d}s+\int_{0}^{t}\widetilde{r}_{s}\,\mathrm{d}B_{s},$
(32)

with $B$ a one-dimensional Brownian motion independent of $W$, and
$\bar{r},\widetilde{r}$ stochastic processes adapted to the natural filtration
of $B$. These dynamics in particular cover the deterministic case
$\widetilde{r}_{t}=0$, so that one can for instance consider time-periodic
fluxes

$\bar{r}_{t}=\sin(2\pi\omega t),\qquad\widetilde{r}_{t}=0,$

or a stochastic process $\mathcal{R}^{r}$ whose ergodic properties are well-
understood, such as an Ornstein–Uhlenbeck process centered at $r$:

$d\mathcal{R}^{r}_{t}=\gamma(r-\mathcal{R}^{r}_{t})\,\mathrm{d}t+s\,\mathrm{d}B_{t}.$

One expects the resulting steady-state to be non-singular with respect to the
invariant measure of the reference dynamics (2), which may be of some use from
a theoretical perspective. The analog of (11) is given by

$\left\\{\begin{aligned}
\mathrm{d}Y^{\mathcal{R}}_{t}&=b(Y^{\mathcal{R}}_{t})\,\mathrm{d}t+\sigma(Y_{t}^{\mathcal{R}})\,\mathrm{d}W_{t}+F(Y_{t}^{\mathcal{R}})\,\mathrm{d}\Lambda^{\mathcal{R}}_{t},\\\
R(Y^{\mathcal{R}}_{t})&=R(Y^{\mathcal{R}}_{0})=\mathcal{R}^{r}_{t}.\end{aligned}\right.$
(33)

Following the same strategy as in the constant response case, one can express
the dynamics in closed form as

$\mathrm{d}Y^{r}_{t}=\overline{P}_{F,\nabla
R}(Y^{\mathcal{R}}_{t})\left[b(Y^{\mathcal{R}}_{t})\,\mathrm{d}t+\sigma(Y^{\mathcal{R}}_{t})\,\mathrm{d}W_{t}\right]+\frac{\bar{r}_{t}\,\mathrm{d}t+\widetilde{r}_{t}\,\mathrm{d}B_{t}}{\nabla
R(Y^{\mathcal{R}}_{t})\cdot F(Y^{\mathcal{R}}_{t})}F(Y_{t}^{\mathcal{R}})\\\
-\frac{1}{2\nabla R(Y^{\mathcal{R}}_{t})\cdot
F(Y^{\mathcal{R}}_{t})}\left(\nabla^{2}R:\left[\frac{F\otimes F}{(\nabla
R\cdot F)^{2}}\widetilde{r}_{t}^{2}+\Pi_{F,\nabla
R,\sigma}\right]\right)(Y^{\mathcal{R}}_{t})F(Y^{\mathcal{R}}_{t})\,\mathrm{d}t,$
(34)

and the bounded-variation contribution to the forcing as

$\lambda^{\mathcal{R}}_{t}=\left[\frac{1}{\nabla R\cdot
F}\left(\bar{r}_{t}-b\cdot\nabla R-\frac{1}{2}\nabla^{2}R:\left[\frac{F\otimes
F}{(\nabla R\cdot F)^{2}}\widetilde{r}_{t}^{2}+\Pi_{F,\nabla
R,\sigma}\right]\right)\right](Y^{\mathcal{R}}_{t}),$

which is still an Itô process adapted to the larger filtration
$\left(\sigma(B_{s},W_{s}:0\leqslant s\leqslant t)\right)_{t\geqslant 0}.$ We
refer to the Appendix A.2 for details of the computations.

##  4 Mobility and shear viscosity computations for Langevin dynamics

So far, we have made very few assumptions about the type of reference
dynamics, driving force or flux. We now turn to presenting a framework in
which the computation of physically meaningful transport coefficients may be
performed, namely that of non-equilibrium (kinetic or underdamped) Langevin
dynamics. We first recall the NEMD method for underdamped Langevin dynamics in
Section 4.1, before specifying the expressions of the non-equilibrium forcings
and response observables relevant for mobility and shear viscosity
computations in Section 4.2. We finally derive the corresponding Norton
dynamics in Section 4.3, and interpret it in terms of a principle of least
constraint.

###  4.1 Standard non-equilibrium Langevin dynamics

The non-equilibrium framework we consider is that of a perturbation of the
Langevin dynamics by an external driving force $F$. The state of the system is
described by a vector $q\in\mathcal{D}$, describing the position of atoms, and
a vector of corresponding momenta $p\in\mathbb{R}^{dN}$, where $d$ is the
physical dimension and $N$ is the number of atoms. The configurational domain
$\mathcal{D}$ is usually considered to be $\mathbb{R}^{dN}$ or
$(L\mathbb{T})^{dN}$ for some box length $L>0$. The evolution of the system is
governed by the Hamiltonian

$H(q,p)=V(q)+\frac{1}{2}p\cdot M^{-1}p,$

where $V:\mathcal{D}\to\mathbb{R}$ is the potential energy function. In
practice, $V$ is determined empirically to give an approximation of the
ground-state energy for the Schrödinger Hamiltonian corresponding to the
Born–Oppenheimer description of the system. The dynamics is given by the
following stochastic differential equation:

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=M^{-1}p_{t}\,\mathrm{d}t,\\\
\mathrm{d}p_{t}&=\left[-\nabla V(q_{t})+\eta F(q_{t})\right]\mathrm{d}t-\gamma
M^{-1}p_{t}\,\mathrm{d}t+\sigma\,\mathrm{d}W_{t}.\end{aligned}\right.$ (35)

The parameter $\beta>0$ is proportional to the inverse temperature, $M$ is a
positive-definite symmetric mass matrix (typically a diagonal matrix), and
$\gamma$ and $\sigma$ are $dN\times dN$ matrices satisfying the fluctuation-
dissipation relation

$\sigma\sigma^{\intercal}=\frac{2\gamma}{\beta}.$

The friction coefficient $\gamma$ is thus a symmetric positive semi-definite
matrix, and $W$ is a standard $dN$-dimensional Brownian motion. Typically,
$\gamma$ is taken to be either a constant or a positive diagonal matrix, so
that one simply writes $\sigma=\sqrt{2\gamma/\beta}$. The parameter $\eta$
once again governs the magnitude of the non-equilibrium perturbation. One can
show that at equilibrium ($\eta=0$), the Boltzmann–Gibbs distribution

$\mu(q,p)\,\mathrm{d}q\,\mathrm{d}p=Z^{-1}\mathrm{e}^{-\beta
H(q,p)}\,\mathrm{d}q\,\mathrm{d}p$ (36)

is an invariant probability measure for (35). In view of the separability of
the Hamiltonian into a configurational and a kinetic part, the Boltzmann–Gibbs
measure can be written in tensor form, as the product of a configurational
measure

$\mathit{m}(q)\,\mathrm{d}q=\frac{1}{Z_{m}}\mathrm{e}^{-\beta
V(q)}\,\mathrm{d}q,$

and of a Gaussian kinetic measure

$\kappa(p)\,\mathrm{d}p=\operatorname{det}\left(\frac{2\pi
M}{\beta}\right)^{-1/2}\mathrm{e}^{-\frac{\beta}{2}p\cdot
M^{-1}p}\,\mathrm{d}p.$

Momenta and positions are thus uncorrelated at equilibrium. Since the
perturbation is generally non-gradient, there is however no way to write out
the expression of the invariant probability measure of the non-equilibrium
dynamics ($\eta\neq 0$). Momenta and positions typically have a correlation of
order $\eta$ under this measure.

###  4.2 Non-equilibrium forcings and fluxes

In this work, we consider transport coefficients corresponding to fluxes which
can be written under the form

$R(q,p)=G(q)\cdot p,$ (37)

where $G:\mathcal{D}\to\mathbb{R}^{dN}$ is a vector field. Such forms of the
response can be understood as measuring correlations between the momenta and
some possibly non-linear feature of the configurational coordinates. In
particular, such responses have zero-mean at equilibrium, owing to the tensor
form of the Boltzmann–Gibbs measure and the fact that the average momentum
under $\kappa$ vanishes. The form (37) is general enough to capture the cases
of mobility and shear viscosity computations for molecular fluids, which we
now proceed to present, and which are our numerical test cases.

#### Mobility computations.

We first describe the NEMD method for the computation of diffusion properties.
Here, we consider a periodic domain $\mathcal{D}=(L\mathbb{T})^{d}$. The NEMD
method is obtained by taking $F$ as a constant vector field in (35), and the
velocity in the direction $F$ as the response, which measures the particle
flux in the direction $F$. We assume that $F$ is normalized as $\|F\|=1$. Thus
the perturbation and response observable are defined respectively by

$F(q)=F\in\mathbb{R}^{3N},\qquad R(q)=F\cdot M^{-1}p.$ (38)

Using the symmetry of $M$, one may rewrite $R(q,p)$ under the form (37). For
practical computations, we consider two cases:

  * $\circ$

_Single drift_ : this corresponds to a perturbation where the force acts on a
single component of the momentum, which can be assumed (by
indistinguishability of the particles) to be the $x$ component of the first
particle:

$F_{\mathrm{S}}=(1,0,\dots)^{\intercal}\in\mathbb{R}^{dN}.$ (39)

  * $\circ$

_Color drift_ (see [14, Chapter 6]): this corresponds to a perturbation in
which the force acts on half of the particles in one direction, and on half of
the particles in the opposite direction, which we choose by convention to be
the $x$ direction:

$F_{\mathrm{C}}=\frac{1}{\sqrt{N}}(\underbrace{1,0,\dots,0}_{\in\,\mathbb{R}^{d}},-1,0,\dots,1,0\dots)^{\intercal}\in\mathbb{R}^{dN}.$
(40)

The corresponding transport coefficients are related to the diffusion
properties of the molecular system. More precisely, the transport coefficient
$\alpha_{F_{\mathrm{S}}}$ for the single drift forcing can be related to the
diffusion coefficient

$D=\underset{T\to\infty}{\lim}\,\frac{1}{2dNT}\sum_{i=1}^{dN}\mathbb{E}_{0}\left[\left(\int_{0}^{T}M^{-1}p_{s,i}\,\mathrm{d}s\right)^{2}\right]$

for an isotropic system via

$\alpha_{F_{\mathrm{S}}}=\beta D,$ (41)

as shown in [38]. When the potential energy function satisfies

$\sum_{i=1}^{N}\frac{\partial}{\partial q_{i,x}}V(q)=0,$ (42)

a condition satisfied for instance for pairwise interactions depending only on
the relative distances (as for (83) below), then the transport coefficient
$\alpha_{F_{\mathrm{C}}}$ for the color drift is related to
$\alpha_{F_{\mathrm{S}}}$ through

$\alpha_{F_{\mathrm{C}}}=\alpha_{F_{\mathrm{S}}}-\frac{2\lfloor
N/2\rfloor}{N(N-1)}\left(\frac{1}{\gamma}-\alpha_{F_{\mathrm{S}}}\right).$
(43)

See Appendix A.3 for the proof of this relation. In particular
$\alpha_{F_{\mathrm{C}}}$ coincides with $\alpha_{F_{\mathrm{S}}}$ in the
thermodynamic limit $N\to\infty$.

The interest of considering these two drift perturbations is to assess whether
the Norton method requires the forcing to act on the bulk of the sytem (color
drift case) to be consistent with the NEMD method. The other extreme is that
of a forcing acting on a single particle (single drift case).

#### Shear viscosity computations.

The Norton framework described in this section also allows for the computation
of the shear viscosity; see [43] for a review of various approaches to compute
this transport coefficient. Instead of considering bulk forcings which require
modifying the periodic boundary conditions (as studied in [10] for instance),
we consider the method introduced in [17], for which the non-gradient force is
periodic. The mathematical properties of this approach in the case of Langevin
dynamics are studied in [24]. In this setting, we consider an anisotropic
three-dimensional configurational domain of the form

$\mathcal{D}=\left(L_{x}\mathbb{T}\times L_{y}\mathbb{T}\times
L_{z}\mathbb{T}\right)^{N}.$ (44)

We further allow for an anisotropic friction coefficient $\gamma$ which is
diagonal and defined by three directional friction coefficients
$\gamma_{x},\gamma_{y},\gamma_{z}\,>0$. A forcing is applied on the momenta in
the longitudinal direction $x$, with an intensity depending on the transverse
configurational coordinate $y$, according to a predefined forcing profile
$f:~{}L_{y}\mathbb{T}~{}\to~{}\mathbb{R}$:

$\forall\,1\leqslant j\leqslant N,\qquad F(q)_{j,x}=f(q_{j,y}),\qquad
F(q)_{j,y}=F(q)_{j,z}=0.$ (45)

In the non-equilibrium steady-state, the system displays an average
longitudinal velocity profile depending only the transverse coordinate and the
magnitude of the non-equilibrium perturbation. More precisely, given an
approximation of the identity
$\left(\psi_{\varepsilon}\right)_{\varepsilon>0}$ on $L_{y}\mathbb{T}$, define

$u_{x}(y)=\underset{\varepsilon\to 0}{\lim}\,\underset{\eta\to
0}{\lim}\,\frac{L_{y}}{\eta
N}\,\mathbb{E}_{\eta}\left[\sum_{j=1}^{N}\left(M^{-1}p\right)_{j,x}\psi_{\varepsilon}(q_{j,y}-y)\right].$
(46)

The term $L_{y}$ in the numerator is motivated by the fact that the average
with respect to the $y$ coordinate of the term inside the limits is the
velocity, up to a factor $1/\eta$:

$\displaystyle\frac{1}{L_{y}}\int_{L_{y}\mathbb{T}}\left(\frac{L_{y}}{\eta
N}\mathbb{E}_{\eta}\left[\sum_{j=1}^{N}\left(M^{-1}p\right)_{j,x}\psi_{\varepsilon}(q_{j,y}-y)\right]\right)\,\mathrm{d}y$
$\displaystyle=\frac{1}{\eta
N}\sum_{j=1}^{N}\mathbb{E}_{\eta}\left[\left(M^{-1}p\right)_{j,x}\right]$
$\displaystyle=\frac{1}{\eta}\mathbb{E}_{\eta}\left[\left(M^{-1}p\right)_{1,x}\right],$

using the indistinguishability of the particles to obtain the last equality.
Thus $u_{x}$ corresponds to a localized linear response of the longitudinal
velocity, which can be estimated in practice from trajectory averages using a
binning procedure. One can then derive the following differential equation
relating $u_{x}$ to the shear viscosity $\nu$:

$-\nu u_{x}^{\prime\prime}(y)+\gamma_{x}\rho u_{x}(y)=\rho f(y),$ (47)

where $\rho=N/(L_{x}L_{y}L_{z})$ is the particle density of the system. Since
this profile is periodic in the transverse coordinate, the shear viscosity can
be related to the Fourier coefficients of $u_{x}$ and $f$, namely

$U_{1}=\frac{1}{L_{y}}\int_{0}^{L_{y}}u_{x}(y)\,\mathrm{e}^{\frac{2\mathrm{i}\pi
y}{L_{y}}}\,\mathrm{d}y,\qquad
F_{1}=\frac{1}{L_{y}}\int_{0}^{L_{y}}f(y)\,\mathrm{e}^{\frac{2\mathrm{i}\pi
y}{L_{y}}}\,\mathrm{d}y,$

through

$\nu=\rho\left(\frac{F_{1}}{U_{1}}-\gamma_{x}\right)\left(\frac{L_{y}}{2\pi}\right)^{2}.$
(48)

For practical purposes, we choose $f$ such that $F_{1}$ is analytically known,
which is for instance the case if $f$ is a sinusoidal profile. Taking the
limit $\varepsilon\to 0$, in (46), the Fourier coefficient $U_{1}$ can be
rewritten as

$U_{1}=\lim_{\eta\to 0}\frac{1}{\eta
N}\mathbb{E}_{\eta}\left[\sum_{j=1}^{N}\left(M^{-1}p\right)_{j,x}\mathrm{exp}\left(\frac{2\mathrm{i}\pi
q_{j,y}}{L_{y}}\right)\right].$ (49)

This expression is precisely the linear response (6) for the following
response observable, which is an empirical Fourier coefficient for the
longitudinal velocity profile:

$R(q,p)=\frac{1}{N}\sum_{j=1}^{N}\left(M^{-1}p\right)_{j,x}\mathrm{exp}\left(\frac{2\mathrm{i}\pi
q_{j,y}}{L_{y}}\right).$ (50)

This response is also of the form (37). One can in turn estimate it using
trajectory averages, as in (7), yielding the estimator

$\widehat{U}_{1,\eta,T}=\frac{1}{\eta
TN}\sum_{j=1}^{N}\int_{0}^{T}\left(M^{-1}p_{t}\right)_{j,x}\exp\left(\frac{2\mathrm{i}\pi
q_{j,y,t}}{L_{y}}\right)\mathrm{d}t.$ (51)

###  4.3 The Norton method for Langevin dynamics

We write more explicitly in this section the Norton dynamics for observables
of the form (37). The general dynamics (11) greatly simplifies in this
situation, and also admits a nice physical interpretation.

#### Closed form of the dynamics.

Note first that $\nabla^{2}_{p}R(q,p)=0$ when $R$ is of the form (37), which
together with the fact that the noise is degenerate (it acts only on the
momenta in (35)) implies that the quadratic covariation term in (21) vanishes.
The forcing $F$, too, acts only on the momenta, so that computing the
projectors (18) for the choices

$A(q,p)=\begin{pmatrix}0\\\ F(q)\end{pmatrix},\qquad B(q,p)=\nabla
R(q,p)=\begin{pmatrix}\nabla G(q)p\\\ G(q)\end{pmatrix},$

yields the following expression for the projector associated with the Norton
dynamics:

$\mathrm{I}_{2dN}-\frac{A\otimes B}{A\cdot
B}(q,p)=\begin{pmatrix}\displaystyle{\mathrm{I}_{dN}}&0\\\
\displaystyle{-\frac{F(q)\otimes\nabla G(q)p}{F(q)\cdot
G(q)}}&\displaystyle{\mathrm{I}_{dN}-\frac{F(q)\otimes G(q)}{F(q)\cdot
G(q)}}\end{pmatrix}.$

Since

$\mathrm{I}_{dN}-\frac{F(q)\otimes G(q)}{F(q)\cdot
G(q)}=\overline{P}_{F,G}(q),$

where we slightly abuse the notation (18), it follows that the Norton dynamics
(21) is given by

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=M^{-1}p_{t}\,\mathrm{d}t,\\\
\mathrm{d}p_{t}&=\overline{P}_{F,G}(q_{t})\left(-\nabla
V(q_{t})\,\mathrm{d}t-\gamma
M^{-1}p_{t}\,\mathrm{d}t+\sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_{t}\right)\\\
&\quad-\frac{\nabla G(q_{t})p_{t}\cdot M^{-1}p_{t}}{F(q_{t})\cdot
G(q_{t})}F(q_{t})\,\mathrm{d}t.\end{aligned}\right.$ (52)

Note that the configurational part of the dynamics is unaffected, which is
consistent with the choice of a non-equilibrium perturbation acting solely on
the momenta. The forcing observable (20) is given by

$\lambda(q,p)=\frac{G(q)\cdot(\nabla V(q)+\gamma M^{-1}p)-\nabla G(q)p\cdot
M^{-1}p}{F(q)\cdot G(q)},$ (53)

and the controllability condition is given by

$(G\cdot F)(q)\neq 0.$ (54)

######  Remark 1.

We have restricted ourselves to a simple Langevin dynamics where the drift
term is the sum of a gradient force and a linear friction term, and the noise
is additive. The derivation can however be extended verbatim to a more general
dynamics of the form

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=M^{-1}p_{t}\,\mathrm{d}t,\\\
\mathrm{d}p_{t}&=b(q_{t},p_{t})\,\mathrm{d}t+\sigma(q_{t},p_{t})\,\mathrm{d}W_{t},\end{aligned}\right.$

for general drift terms $b$ and diffusion matrices $\sigma$. For these kinetic
dynamics, the associated Norton dynamics is given by

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=M^{-1}p_{t}\,\mathrm{d}t,\\\
\mathrm{d}p_{t}&=\overline{P}_{F,G}(q_{t})\left[b(q_{t},p_{t})\,\mathrm{d}t+\sigma(q_{t},p_{t})\,\mathrm{d}W_{t}\right]-\frac{\nabla
G(q_{t})p_{t}\cdot M^{-1}p_{t}}{F(q_{t})\cdot
G(q_{t})}F(q_{t})\,\mathrm{d}t.\end{aligned}\right.$

#### Physical interpretation.

In the particular case of Langevin dynamics, the Norton method has the
physical interpretation that it obeys a version of Gauss’s principle of least
constraint. The connection between Gauss’s principle and non-equilibrium
thermodynamics was already pointed out in [13], see also the discussion in
[14, Section 5.2]. The principle of least constraints is a statement of
classical dynamics equivalent to D’Alembert’s principle, stating that the
force applied to a system subject to a set of holonomic or non-holonomic
constraints minimizes at every point in time the Euclidean distance to the
force of the same system free from any constraints. Although we describe this
interpretation in a deterministic setting, we stress that it remains valid, at
least on a formal level, upon considering a stochastic version of the
dynamics. More precisely, we assume that the dynamics for the unconstrained
system can be written under the form

$\left\\{\begin{aligned} \dot{q}&=M^{-1}p,\\\
\dot{p}&=f_{\mathrm{ref}}(q,p),\end{aligned}\right.$ (55)

and that the constraint is of the form $R(q)=r$ (holonomic case) or $R(q,p)=r$
(non-holonomic case). We can then write the dynamics of the constrained system
as

$\left\\{\begin{aligned} \dot{q}&=M^{-1}p,\\\
\dot{p}&=f_{\mathrm{cons}}(q,p),\end{aligned}\right.$ (56)

where $f_{\mathrm{cons}}(q,p)$ is the force on the constrained system obeying
Gauss’s principle, which dictates that $f_{\mathrm{cons}}(q,p)$ is the
orthogonal projection of $f_{\mathrm{ref}}(q)$ onto the affine hyperplane
$\mathcal{H}_{q,p}$ of admissible forces. This hyperplane can be determined by
differentiating the constraint in time and setting it to zero. To obtain a
constraint on the time derivative of the momenta $\dot{p}$, this
differentiation has to be performed once in the non-holonomic case and twice
in the holonomic case. A simple computation shows that the hyperplane of
admissible forces is thus given in the holonomic case by

$\mathcal{H}_{q,p}=\left\\{\xi\in\mathbb{R}^{dN}\,\middle|\,\xi\cdot
M^{-1}\nabla R(q)+\left(M^{-1}p\right)^{\otimes 2}:\nabla^{2}R(q)=0\right\\},$
(57)

and in the non-holonomic case by

$\mathcal{H}_{q,p}=\left\\{\xi\in\mathbb{R}^{dN}\,\middle|\,\xi\cdot\nabla_{p}R(q,p)+\left(M^{-1}p\right)\cdot\nabla_{q}R(q,p)=0\right\\}.$
(58)

In particular, $f_{\mathrm{cons}}(q,p)-f_{\mathrm{ref}}(q,p)$ is orthogonal to
$\nabla_{p}R(q,p)$ in the non-holonomic case, and to $M^{-1}\nabla R(q)$ in
the holonomic case. The use of holonomic constraints is widespread in MD,
where they are used to simulate systems with molecular constraints (such as
fixed bond lengths or bond angles). See [29, Chapter 4] for a general
introduction, and [33] for a detailed study of the underdamped Langevin case.

The Norton dynamics (52) for a response function $R$ of the form (37) also
satisfies a version of Gauss’s principle of least constraint with respect to a
non-Euclidean metric for which $F(q)$ is everywhere orthogonal to
$G(q)^{\perp}$. More precisely, this metric is induced by the configuration-
dependent norm

$\|\xi\|^{2}_{q}=\left(\xi\cdot F(q)\right)^{2}+\left\|\xi-\frac{\xi\cdot
G(q)}{\|G(q)\|^{2}}G(q)\right\|^{2},$ (59)

where $\|\cdot\|$ denotes the usual Euclidean norm. This metric is defined and
non-degenerate if $F$ and $G$ are both non-zero, which is implied by the
controllability condition (22). As we show below, the force on the constrained
system is, in this metric, obtained as

$\underset{\xi\in\mathcal{H}_{q,p}}{\operatorname*{argmin}}\left\|\xi-
f_{\mathrm{ref}}(q,p)\right\|_{q}=\overline{P}_{F,G}(q)f_{\mathrm{ref}}(q,p)-\frac{\nabla
G(q)p\cdot M^{-1}p}{F(q)\cdot G(q)}F(q),$ (60)

where $\mathcal{H}_{q,p}$ is the hyperplane defined in (58) associated with
the non-holonomic constraint $R(q,p)=G(q)\cdot p=r$. This is precisely the
Norton force corresponding to the unconstrained system (55).

To prove (60), we proceed as follows. Note first that $\mathcal{H}_{q,p}$ in
(58) is an affine translate of $\nabla_{p}R(q,p)^{\perp}=G(q)^{\perp}$, the
Euclidean orthogonal to $G(q)$, so that its normal direction with respect to
the scalar product induced by $\|\cdot\|_{q}$ is $F(q)$. Furthermore, the norm
$\|\cdot\|_{q}$ is constructed precisely so that the projector
$\overline{P}_{F,G}(q)$, whose action is depicted in Figure 1, is an
orthogonal projector onto $G(q)^{\perp}$, for the scalar product associated
with $\|\cdot\|_{q}$. This implies that the minimizer in (60) is of the form

$\xi^{*}(q,p)=f_{\mathrm{ref}}(q,p)-\alpha F(q)\in\mathcal{H}_{q,p}.$

The value of $\alpha$ is determined by the condition

$\left(f_{\mathrm{ref}}(q,p)-\alpha F(q)\right)\cdot G(q)+\nabla G(q)p\cdot
M^{-1}p=0,$

hence

$\alpha=\frac{f_{\mathrm{ref}}(q)\cdot G(q)+\nabla G(q)p\cdot
M^{-1}p}{F(q)\cdot G(q)},$

and finally

$\xi^{*}(q,p)=f_{\mathrm{ref}}(q)-\alpha
F(q)=\overline{P}_{F,G}(q)f_{\mathrm{ref}}(q)-\frac{\nabla G(q)p\cdot
M^{-1}p}{F(q)\cdot G(q)}F(q).$

Loosely speaking, the Norton dynamics is the least non-equilibrium like of all
dynamics on the constant response manifold, if one measures similarity in
terms of the force considered in the $\|\cdot\|_{q}$ metric.

##  5 Numerical discretizations of Norton dynamics

We describe in this section a discretization of the Norton dynamics, first for
the general dynamics (11) in Section 5.1, before specializing it to the
setting of Langevin dynamics. The numerical schemes are based for Langevin
dynamics on splitting schemes, inspired by the ones typically used in the NEMD
setting, but have the additional property of exactly preserving the flux
throughout the numerical trajectory. For completeness, we describe splitting
schemes in the case of standard NEMD Langevin dynamics in Section 5.2, and
then in the Norton setting in Section 5.3. We finally describe in Section 5.4
a way to estimate the bounded variation part of the forcing (20) from the
numerical trajectories, in order to approximate (25).

###  5.1 Numerical schemes for general Norton dynamics

We discuss here the simulation of the Norton dynamics (11), which requires a
discretization in time. Formally, given a time step $\Delta t>0$, a numerical
scheme for the equilibrium dynamics (2) is defined by a map $\Phi_{\Delta
t}:\mathcal{X}\times\mathbb{R}^{m}\to\mathcal{X}$ which, iterated with
independent and identically distributed (i.i.d.) standard $m$-dimensional
Gaussian variables $(\mathcal{G}_{n})_{n\geqslant 0}$ yields a discretization
of the dynamics:

$X^{n+1}=\Phi_{\Delta t}(X^{n},\mathcal{G}^{n}),$ (61)

i.e. $X^{n}$ is an approximation of $X_{n\Delta t}$. For example, the
Euler–Maruyama scheme for the Thévenin dynamics (2) corresponds to

$\Phi^{\mathrm{EM}}_{\Delta t}(x,g)=x+\Delta t\left[b(x)+\eta
F(x)\right]+\sqrt{\Delta t}\,\sigma(x)g.$ (62)

In principle, it would be possible to consider discretizations directly based
on the autonomous form of the Norton dynamics (21), and average the forcing
observable $\lambda$ along the so-obtained numerical trajectories. However, a
key property we require from a numerical scheme is that the constant flux
manifold should be preserved under the discrete dynamics. Standard
discretizations of (21) usually do not satisfy such preservation properties.
Moreover, since the autonomous form of the dynamics involves second-order
derivatives of the response, which may be cumbersome or expensive to compute,
it is typically more convenient numerically to take another approach, obtained
by enforcing the constraint via a Lagrange multiplier. Given a stochastic
scheme $\Phi_{\Delta t}$ for the reference dynamics, we can consider the
following discretization of the Norton dynamics:

$\left\\{\begin{aligned} \widetilde{X}^{n+1}&=\Phi_{\Delta t}(X^{n},G^{n}),\\\
X^{n+1}&=\widetilde{X}^{n+1}+\Delta t\Lambda^{n,*}F(X^{n}),\\\
R(X^{n+1})&=r.\end{aligned}\right.$ (63)

Note that we chose here, somewhat arbitrarily, to perform the projection with
respect to $F(X^{n})$. Other possible choices include $F(\widetilde{X}^{n+1})$
or $F(X^{n+1})$, the latter choice generally yielding an implicit scheme. In
any case, finding $\Lambda^{n,*}$ requires solving

$R\left(\widetilde{X}^{n+1}+\Delta t\Lambda^{n,*}F(X^{n})\right)=r,$

which is typically a non-linear equation, for which the appropriate numerical
strategy depends on the situation at hand. A typical choice is to resort to
Newton’s method or a fixed-point iteration, as done when enforcing holonomic
constraints in MD (see [20, Section VII.1] and [31, Chapter 7]), as well as
[40, 2] for pioneering works motivated by applications in MD.

Note that $\Lambda^{n,*}$ corresponds to an approximation of the full forcing
increment $\Lambda_{(n+1)\Delta t}-\Lambda_{n\Delta t}$, and in particular
incorporates the martingale increment. This martingale part should be removed
when estimating $\mathbb{E}_{r}^{*}[\lambda]$ in order to reduce the variance
of the estimator under consideration. This can be done at first order in
$\sqrt{\Delta t}$ by using a control variate derived from its analytic
expression (16), and the equality in law $\sqrt{\Delta
t}\,G^{n}\overset{\mathrm{law}}{=}W_{(n+1)\Delta t}-W_{n\Delta t}$. This leads
to the following estimator of $\lambda\Delta t$:

$\Lambda^{n}=\Lambda^{n,*}+\sqrt{{\Delta t}}\frac{\nabla
R(X^{n})\cdot\sigma(X^{n})G^{n}}{\nabla R(X^{n})\cdot F(X^{n})}.$ (64)

The efficiency of such variance reduction methods has already been
demonstrated for the estimation of mean Lagrange multipliers associated with
holonomic constraints, for example in the context of free energy computations
(see [8] and [32, Remark 3.33]). An estimation of the ensemble average
$\mathbb{E}_{r}^{*}[\lambda]$ can be obtained by considering trajectory
averages

$\widehat{\lambda}_{N_{\mathrm{iter}}}=\frac{1}{N_{\mathrm{iter}}}\sum_{n=1}^{N_{\mathrm{iter}}}\lambda^{n},\qquad\lambda^{n}=\frac{1}{{\Delta
t}}\Lambda^{n}.$ (65)

As hinted at above, this strategy has the clear advantage that one does not
have to compute $\lambda(X^{n})$ along the trajectory (which may be tedious to
implement in practice), since consistent approximations thereof appear as
natural byproducts of the integration procedure.

###  5.2 Splitting schemes for (non-)equilibrium Langevin dynamics

In the particular case where the reference dynamics is the Langevin dynamics
(35), one can resort to a class of discretization strategies based on operator
splittings of the infinitesimal generator, which we briefly recall for
completeness in this section (see [27, 30], and [29, Chapter 7] for more
details). By carefully choosing the order of the operators appearing in this
splitting in the equilibrium setting, one can devise highly stable numerical
schemes, which have a lower bias when sampling configurational averages in the
overdamped limit $\gamma\to\infty$ compared with other choices of the
splitting order (see [28]), while also showing desirable energy-conservation
properties in the Hamiltonian limit $\gamma\to 0$. Such schemes have become
the standard choice to integrate kinetic Langevin dynamics in MD applications.
The general procedure is as follows.

The generator of the (non-)equilibrium Langevin dynamics can be written as the
sum of four terms,

$\mathcal{L}_{\eta}=\mathcal{L}^{\mathrm{A}}+\mathcal{L}^{\mathrm{B}}+\gamma\mathcal{L}^{\mathrm{O}}+\eta\widetilde{\mathcal{L}},$

with

$\left\\{\begin{aligned} \mathcal{L}^{\mathrm{A}}&=M^{-1}p\cdot\nabla_{q},\\\
\mathcal{L}^{\mathrm{B}}&=-\nabla V\cdot\nabla_{p},\\\
\mathcal{L}^{\mathrm{O}}&=-M^{-1}p\cdot\nabla_{p}+\frac{1}{\beta}\Delta_{p},\\\
\widetilde{\mathcal{L}}&=F\cdot\nabla_{p}.\end{aligned}\right.$ (66)

The operators $\mathcal{L}^{\mathrm{A}}$,
$\mathcal{L}^{\mathrm{B}}+\eta\widetilde{\mathcal{L}}$ and
$\gamma\mathcal{L}^{\mathrm{O}}$ can be viewed as infinitesimal generators in
their own right, which correspond respectively to the three following
elementary dynamics. The dynamics generated by $\mathcal{L}^{\mathrm{A}}$ is
given by

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=M^{-1}p_{t}\,\mathrm{d}t,\\\
\mathrm{d}p_{t}&=0,\end{aligned}\right.$ (67)

while $\mathcal{L}^{\mathrm{B}}+\eta\widetilde{\mathcal{L}}$ gives rise to

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=0,\\\ \mathrm{d}p_{t}&=\left(-\nabla
V(q_{t})+\eta F(q_{t})\right)\,\mathrm{d}t.\end{aligned}\right.$ (68)

Finally, $\mathcal{L}^{\mathrm{O}}$ generates the following Ornstein–Uhlenbeck
process on the momenta:

$\left\\{\begin{aligned} \mathrm{d}q_{t}&=0,\\\ \mathrm{d}p_{t}&=-\gamma
M^{-1}p_{t}\,\mathrm{d}t+\sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_{t}.\end{aligned}\right.$
(69)

The three elementary dynamics (67), (68), (69) are analytically integrable.
Splitting schemes for the Thévenin dynamics are obtained by composing the
evolution operators $\mathrm{e}^{{\Delta t}\mathcal{L}^{\mathrm{A}}}$,
$\mathrm{e}^{{\Delta
t}\left(\mathcal{L}^{\mathrm{B}}+\eta\widetilde{\mathcal{L}}\right)}$ and
$\mathrm{e}^{{\Delta t}\mathcal{L}^{\mathrm{O}}}$ corresponding to each of
these elementary dynamics. For instance, the evolution operator for the
celebrated $\mathrm{BAOAB}$ method corresponds to

$\mathrm{e}^{\frac{{\Delta
t}}{2}\left(\mathcal{L}^{\mathrm{B}}+\eta\widetilde{\mathcal{L}}\right)}\mathrm{e}^{\frac{{\Delta
t}}{2}\mathcal{L}^{\mathrm{A}}}\mathrm{e}^{{\Delta
t}\gamma\mathcal{L}^{\mathrm{O}}}\mathrm{e}^{\frac{{\Delta
t}}{2}\mathcal{L}^{\mathrm{A}}}\mathrm{e}^{\frac{{\Delta
t}}{2}\left(\mathcal{L}^{\mathrm{B}}+\eta\widetilde{\mathcal{L}}\right)}.$

These schemes can be formally justified, and rigorously analyzed with the
Baker–Campbell–Hausdorff formula, yielding error estimates à la Talay–Tubaro
on the invariant measure in the linear response regime where $|\eta|$ is
small, as well as on estimators of transport coefficients, see [30].

###  5.3 Splitting schemes for Langevin–Norton dynamics

For the Langevin dynamics (52) in the Norton setting, a strategy similar to
the one in Section 5.2 can be used. The generator of the Norton dynamics can
be split as the sum of three operators, each of which is the generator of an
elementary dynamics, still denoted by $\mathrm{A}$, $\mathrm{B}$ and
$\mathrm{O}$. In fact, these are the Norton counterparts of the elementary
dynamics (67), (68) and (69), which therefore preserve the flux by
construction. This splitting strategy is still motivated by the fact that each
elementary dynamics individually preserves the flux observable (37). However,
in contrast to the NEMD case, the flow of the $\mathrm{A}$-dynamics is
generally not known in analytical form, and one has to resort to a numerical
scheme to approximate it.

In order to make the numerical scheme precise, the first step is to split the
generator as

$\mathfrak{L}=\mathfrak{L}^{\mathrm{A}}+\mathfrak{L}^{\mathrm{B}}+\gamma\mathfrak{L}^{\mathrm{O}},$
(70)

with

$\left\\{\begin{aligned}
\mathfrak{L}^{\mathrm{A}}&=M^{-1}p\cdot\nabla_{q}-\frac{\left(\nabla
G\right)p\cdot M^{-1}p}{F\cdot G}F\cdot\nabla_{p},\\\
\mathfrak{L}^{\mathrm{B}}&=-\overline{P}_{F,G}\nabla V\cdot\nabla_{p},\\\
\mathfrak{L}^{\mathrm{O}}&=-\overline{P}_{F,G}M^{-1}p\cdot\nabla_{p}+\frac{1}{\beta}\overline{P}_{F,G}\overline{P}_{F,G}^{\intercal}:\nabla_{p}^{2}.\end{aligned}\right.$
(71)

This decomposition is motivated by the flux conservation properties
$\mathfrak{L}^{\mathrm{B}}R=\mathfrak{L}^{\mathrm{O}}R=0$, which immediately
imply that $\mathfrak{L}^{\mathrm{A}}R=0$, since the overall dynamics with
generator $\mathfrak{L}$ conserves the non-equilibrium response by
construction.

#### Decomposition into elementary dynamics.

Analogously to the NEMD case, we interpret the various operators on the right-
hand side of (70) as the generators of elementary dynamics, and use the same
terminology as in Section 5.2. In fact, these elementary dynamics are
precisely the Norton counterparts (21) to each of the NEMD elementary
dynamics, so that they indeed individually preserve the response. They read

$\mathrm{A}$ dynamics:  $\displaystyle\left\\{\begin{aligned}
\mathrm{d}q_{t}&=M^{-1}p_{t}\,\mathrm{d}t,\\\ \mathrm{d}p_{t}&=-\frac{\nabla
G(q_{t})p_{t}\cdot M^{-1}p_{t}}{F(q_{t})\cdot
G(q_{t})}F(q_{t})\,\mathrm{d}t.\end{aligned}\right.$ (72) $\mathrm{B}$
dynamics:  $\displaystyle\left\\{\begin{aligned} \mathrm{d}q_{t}&=0,\\\
\mathrm{d}p_{t}&=-\overline{P}_{F,G}(q_{t})\nabla
V(q_{t})\,\mathrm{d}t.\end{aligned}\right.$ $\mathrm{O}$ dynamics:
$\displaystyle\left\\{\begin{aligned} \mathrm{d}q_{t}&=0,\\\
\mathrm{d}p_{t}&=\overline{P}_{F,G}(q_{t})\left(-\gamma
M^{-1}p_{t}\,\mathrm{d}t+\sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_{t}\right).\end{aligned}\right.$

Note that, in contrast to the NEMD setting, the $\mathrm{A}$-dynamics is not
analytically solvable, so that one has to resort to a numerical approximation
similar to those discussed in Section 5.1 to integrate it (see Section 5.4).
However, the Norton counterparts of the $\mathrm{B}$ and $\mathrm{O}$ dynamics
remain analytically integrable, as we discuss below.

#### Analytic integration of the $\mathrm{B}$ and $\mathrm{O}$ dynamics.

The $\mathrm{B}$-dynamics is a time-linear evolution in the momenta, whose
solution is given by

$(q_{t},p_{t})=\left(q_{0},p_{0}-t\overline{P}_{F,G}\nabla
V(q_{0})\right)=\left(q_{0},P_{F,G}(q_{0})p_{0}+\overline{P}_{F,G}(q_{0})\left[p_{0}-t\nabla
V(q_{0})\right]\right).$

In view of the equality $G(q)\cdot p=G(q)\cdot{P}_{F,G}(q)p$, and since
$\overline{P}_{F,G}$ is a projector onto $G^{\perp}$, it is immediate that the
$\mathrm{B}$-dynamics preserves the response flux. One could also simply
notice that it is the Norton counterpart to the NEMD $\mathrm{B}$-dynamics
(68).

The $\mathrm{O}$-dynamics is a projected Ornstein–Uhlenbeck process. To
analytically integrate this dynamics, we assume that $\overline{P}_{F,G}$ and
$\gamma M^{-1}$ commute. If this is not the case, the analytical integration
below should be replaced by an appropriate numerical scheme, such as a
midpoint rule. Using standard arguments of Itô calculus, and the fact that
$\overline{P}_{F,G}$ is a projector, it is straightforward to check that the
solution is given for all $t\geqslant 0$ by $q_{t}=q_{0}$ and

$p_{t}={P}_{F,G}(q_{0})p_{0}+\overline{P}_{F,G}(q_{0})\left(\mathrm{e}^{-t\gamma
M^{-1}}p_{0}+\int_{0}^{t}\mathrm{e}^{-\gamma
M^{-1}(t-s)}\sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_{s}\right).$

This is a Gaussian process, which has the following alternative representation
in distribution:

$p_{t}={P}_{F,G}(q_{0})p_{0}+\overline{P}_{F,G}(q_{0})\left(\mathrm{e}^{-t\gamma
M^{-1}}p_{0}+\sqrt{\frac{1-\mathrm{e}^{-2t\gamma
M^{-1}}}{\beta}M}\mathcal{G}\right),$ (73)

where $\mathcal{G}$ is a standard $dN$-dimensional Gaussian. By the same
arguments as for the $\mathrm{B}$-dynamics, the Norton $\mathrm{O}$-dynamics
preserves the response flux, and is likewise the Norton counterpart to the
NEMD $\mathrm{O}$-dynamics (69).

###  5.4 Estimation of the average forcing

By composing the evolution operators of elementary Norton dynamics, one
obtains a natural splitting approximation of the evolution operator for the
forcing process of the full Norton dynamics. The numerical translation of this
observation is that one can estimate trajectory averages of $\lambda$ directly
from examining individual integration steps of the splitting scheme. We
describe in this section the general procedure to this end.

For a fixed time step $\Delta t>0$ and response magnitude $r\in\mathbb{R}$, we
define three discrete flows acting on the augmented state
$(q,p,\ell)\in\mathcal{D}\times\mathbb{R}^{dN}\times\mathbb{R}$. The role of
the auxiliary variable $\ell$ is to accumulate the bounded-variation component
of Lagrange multipliers enforcing the constant-flux constraint during the
integration step. It is thus initialized at $0$ at the start of each
integration step, accumulating the contributions of each part of the
splitting.

#### Discrete integration of the $\mathrm{A}$ dynamics.

The discrete flow associated with the $\mathrm{A}$-step is given by

${\Phi_{\Delta t,r}^{\mathrm{A}}\left(q,p,\ell\right)=\left(q+\Delta
tM^{-1}p,p+\xi_{\Delta t,r}^{\mathrm{A}}(q,p)F\left(q+\Delta
tM^{-1}p\right),\ell+\xi^{\mathrm{A}}_{\Delta t,r}(q,p)\right)},$ (74)

where $\xi_{\Delta t,r}^{\mathrm{A}}(q,p)$ is the Lagrange multiplier
characterized by the constant-flux condition

$R\left(\Phi_{\Delta_{t},r}^{\mathrm{A}}(q,p,\ell)\right)=r,$

with some abuse of notation in the arguments of $R$. When the response is of
the form $R(q,p)=G(q)\cdot p$, we can explicitly solve for the Lagrange
multiplier, as

$\xi_{\Delta t,r}^{\mathrm{A}}(q,p)=\frac{r-G\left(q+\Delta
tM^{-1}p\right)\cdot p}{F\left(q+\Delta tM^{-1}p\right)\cdot G\left(q+\Delta
tM^{-1}p\right)}.$ (75)

The Norton $\mathrm{A}$-step is therefore equivalent to first evolving the
system’s coordinate variables according to the corresponding NEMD
$\mathrm{A}$-step before correcting the momenta, by applying a constraining
force in the direction of $F(q+\Delta tM^{-1}p)$, the forcing evaluated at the
updated position. As mentioned after (63), there is some freedom in the choice
of point at which to evaluate the forcing when projecting onto the submanifold
$\Sigma_{r}=R^{-1}\\{r\\}$. The motivation for choosing $F(q+\Delta tM^{-1}p)$
is that the two functions appearing in the scalar product $F\cdot G$ in the
denominator of the left-hand side of (75) are evaluated at the same point.

#### Discrete integration of the $\mathrm{B}$ dynamics.

The discrete flow associated with the Norton $\mathrm{B}$-step is given by

$\Phi_{\Delta t,r}^{\mathrm{B}}(q,p,\ell)=\left(q,p-\Delta t\nabla
V(q)+\xi_{\Delta t,r}^{\mathrm{B}}(q,p)F(q),\ell+\xi_{\Delta
t,r}^{\mathrm{B}}(q,p)\right),$ (76)

where $\xi_{\Delta t,r}^{\mathrm{B}}$ is again a Lagrange multiplier enforcing
the constraint, given for a response $G(q)\cdot p$ by

$\xi_{\Delta t,r}^{\mathrm{B}}(q,p)=\frac{r-G\left(q\right)\cdot\left(p-\Delta
t\nabla V(q)\right)}{F\left(q\right)\cdot G\left(q\right)}.$ (77)

This again corresponds to a step of the NEMD $\mathrm{B}$-dynamics, re-
projected onto $\Sigma_{r}$ in the direction $F(q)$. This coincides in fact
with the analytic integration of the elementary $\mathrm{B}$-dynamics (72)
when $G(q)\cdot p=r$.

#### Discrete integration of the $\mathrm{O}$ dynamics.

The flow associated with the Norton $\mathrm{O}$-step is stochastic, hence we
formulate it conditionally on a standard $dN$-dimensional Gaussian
$\mathcal{G}$ increment. In fact, we show below that a convenient update is

${\Phi_{\Delta
t,r}^{\mathrm{O}}\left(q,p,\ell\middle|\,\mathcal{G}\right)=\left(q,\alpha_{\Delta
t}p+\sigma_{\Delta t}\,\sqrt{\Delta t}\,\mathcal{G}+\xi_{\Delta
t,r}^{\mathrm{O}}(q,p,\mathcal{G})F(q),\ell+\frac{r(1-\alpha_{\Delta
t})}{F(q)\cdot G(q)}\right)},$ (78)

where $\alpha_{{\Delta t}}$ and $\sigma_{{\Delta t}}$ are given by

$\alpha_{\Delta t}=\mathrm{e}^{-\gamma M^{-1}{\Delta
t}},\,\qquad\sigma_{\Delta t}=\sqrt{\frac{1-\alpha_{{\Delta
t}}^{2}}{\beta\Delta t}M}.$

To motivate (78), we start by deriving the expression of the Lagrange
multiplier $\xi_{\Delta t,r}^{\mathrm{O}}$, which can be solved for our
particular form of response (37) as

$\xi_{{\Delta t},r}^{\mathrm{O}}(q,p,\mathcal{G})=\frac{r-\alpha_{{\Delta
t}}G(q)\cdot p-\sigma_{{\Delta t}}\sqrt{{\Delta
t}}\,G(q)\cdot{\mathcal{G}}}{F(q)\cdot G(q)}.$ (79)

Again, this corresponds to integrating the NEMD $\mathrm{O}$-dynamics over one
time step before correcting the momenta in the direction $F(q)$. Since the
contribution of $\mathcal{G}$ to this $\xi_{\Delta t,r}^{\mathrm{O}}$ is a
centered Gaussian, we can remove it entirely to only accumulate the non-
martingale component in $\ell$. This is exactly equivalent to the variance
reduction technique discussed in (64). Using $G(q)\cdot p=r$, one arrives at
the expression given in (78) for the action of $\Phi_{{\Delta
t},r}^{\mathrm{O}}$ on $\ell$. Once again, $\Phi_{{\Delta
t},r}^{\mathrm{O}}(q,p,\ell)_{q,p}$ corresponds to the analytic flow of the
Norton $\mathrm{O}$-dynamics (72) over one time step $\Delta t$, with
deterministic initial condition $(q,p)$.

#### Construction of the splitting scheme.

The steps $\mathrm{A,B,O}$ can be composed according to the order of the
splitting, and the bounded-variation increment of the constraining process
over one time step can be estimated from the final increment in $\ell$. For
example, conditionally on $\mathcal{G}^{n+1}$, a standard $dN$ dimensional
Gaussian, the update rule for the Norton BAOAB scheme is given by

${(q^{n+1},p^{n+1},\ell^{n+1})=\Phi_{{\Delta
t}/2,r}^{\mathrm{B}}\circ\Phi_{{\Delta t}/2,r}^{\mathrm{A}}\circ\Phi_{{\Delta
t},r}^{\mathrm{O}}(\cdot\,|\,\mathcal{G}^{n+1})\circ\Phi_{{\Delta
t}/2,r}^{\mathrm{A}}\circ\Phi_{{\Delta t}/2,r}^{\mathrm{B}}(q^{n},p^{n},0)}.$
(80)

Note that we start from $\ell^{n}=0$, and then accumulate the various
increments of the constraining process to obtain $\ell^{n+1}$. For
concreteness, let us explicitly write the update rules for the numerical
scheme (80) and a response observable $R(q,p)=G(q)\cdot p$:

$\left\\{\begin{aligned} p^{n+\frac{1}{5}}&=p^{n}-\frac{\Delta t}{2}\nabla
V(q^{n})+\xi_{\Delta t/2,r}^{\mathrm{B}}(q^{n},p^{n})F(q^{n}),\\\
\ell^{n+\frac{1}{5}}&=\xi_{\Delta t/2,r}^{\mathrm{B}}(q^{n},p^{n}),\\\
q^{n+\frac{2}{5}}&=q^{n}+\frac{\Delta t}{2}M^{-1}p^{n+\frac{1}{5}},\\\
p^{n+\frac{2}{5}}&=p^{n+\frac{1}{5}}+\xi_{\Delta
t/2,r}^{\mathrm{A}}(q^{n},p^{n+\frac{1}{5}})F(q^{n+\frac{2}{5}}),\\\
\ell^{n+\frac{2}{5}}&=\ell^{n+\frac{1}{5}}+\xi_{\Delta
t/2,r}^{\mathrm{A}}(q^{n},p^{n+\frac{1}{5}}),\\\
{p}^{n+\frac{3}{5}}&=\alpha_{\Delta t}p^{n+\frac{2}{5}}+\sigma_{\Delta
t}\mathcal{G}^{n+1}+\xi_{\Delta
t,r}^{\mathrm{O}}(q^{n+\frac{2}{5}},p^{n+\frac{2}{5}},\mathcal{G}^{n+1})F(q^{n+\frac{2}{5}}),\\\
\ell^{n+\frac{3}{5}}&=\ell^{n+\frac{2}{5}}+\frac{r(1-\alpha_{\Delta
t})}{F(q^{n+\frac{2}{5}})\cdot G(q^{n+\frac{2}{5}})},\\\
q^{n+1}&=q^{n+\frac{2}{5}}+\frac{\Delta t}{2}M^{-1}p^{n+\frac{3}{5}},\\\
p^{n+\frac{4}{5}}&=p^{n+\frac{3}{5}}+\xi_{\Delta
t/2,r}^{\mathrm{A}}(q^{n+\frac{2}{5}},p^{n+\frac{3}{5}})F(q^{n+1}),\\\
\ell^{n+\frac{4}{5}}&=\ell^{n+\frac{3}{5}}+\xi_{\Delta
t/2,r}^{\mathrm{A}}(q^{n+\frac{2}{5}},p^{n+\frac{3}{5}}),\\\
{p}^{n+1}&=p^{n+\frac{4}{5}}-\frac{\Delta t}{2}\nabla V(q^{n+1})+\xi_{\Delta
t/2,r}^{\mathrm{B}}(q^{n+1},{p}^{n+\frac{4}{5}})F(q^{n+1}),\\\
\ell^{n+1}&=\ell^{n+\frac{4}{5}}+\xi_{\Delta
t/2,r}^{\mathrm{B}}(q^{n+1},{p}^{n+\frac{4}{5}}).\end{aligned}\right.$ (81)

The average value of $\lambda$ over the corresponding time step can then be
estimated via

$\lambda^{n+1}={\Delta t}^{-1}\ell^{n+1}.$ (82)

Since all the substeps are by construction preserving the value of the
response function, splitting schemes based on (70) define discrete-time Markov
chains on $\Sigma_{r}$. Fixing a splitting and a time step $\Delta t$, one can
hope that the Markov chain admits a unique invariant probability measure
$\pi^{*}_{r,{\Delta t}}$ close to the invariant measure of the continuous time
Norton dynamics for $\Delta t>0$ sufficiently small. Assessing the quality of
these schemes (measured in terms of weak or strong error estimates, or errors
on the invariant measure) as a function of the time step ${\Delta t}$, the
magnitude of the perturbation $r$ and the ordering of the splitting is an open
question.

######  Remark 2.

Let us emphasize that, when the form of the response function does not allow
for analytical expressions of the Lagrange multipliers, one can still apply a
splitting procedure similar to (81), upon replacing the explicit expressions
of $\xi_{\Delta t,r}^{X}$ (for $X\in\\{\mathrm{A},\mathrm{B},\mathrm{O}\\}$)
by the numerical solution of some nonlinear equation determining these
quantities. It is still possible to cancel at dominant order the martingale
increment by subtracting the first order approximation in $\sqrt{\Delta t}$ of
the latter quantity.

##  6 Numerical results

We present in this section the results of various numerical simulations for
the Lennard–Jones fluid. The numerical experiments we perform have several
aims:

  * $\circ$

The first one is to verify that, at least in our setting, the Norton dynamics
is a viable approach for the computation of transport coefficients, in the
sense that Norton estimators of the linear response coincide with those
obtained from NEMD computations. Of course, this cannot be expected for low
dimensional systems. Consider for example the case of kinetic Langevin
dynamics for a single one-dimensional particle on the unit torus $\mathbb{T}$,
with $M=F=1$, and $R(q,p)=p$. The invariant measure for the Norton dynamics is
then easily seen to be the product of the uniform measure on $\mathbb{T}$ with
the Dirac measure $\delta_{r}$ on momentum space, from which it follows in
view of (53) that

$\mathbb{E}_{r}^{*}[\lambda]=\int_{0}^{1}V^{\prime}(q)\,\mathrm{d}q+\gamma
r=\gamma r,$

since $V$ is periodic, whence the Norton transport coefficient is
$\gamma^{-1}$, which differs in general from its NEMD counterpart.

  * $\circ$

Once the validity of the Norton approach has been established, the second aim
is to assess the numerical efficiency of the Norton method, relatively to the
standard NEMD method. A crucial point is to determine whether the Norton
approach leads to a reduction in the asymptotic variance (26) compared to (10)
for estimators of the transport coefficient, as this quantity determines the
magnitude of the statistical error, and therefore the simulation time needed
to attain a given level of accuracy.

In Section 6.1, we present the system and notation used throughout all
numerical experiments. In Section 6.2, we show instances of excellent
agreement in the linear response between Norton and NEMD dynamics. We also
exhibit a case in which the response profiles do not coincide. We further show
that the agreement (when it holds) extends far into the non-linear response
regime. In Section 6.3, we investigate the properties of the Norton and NEMD
systems in the thermodynamic limit, showing that the linear responses coincide
far before convergence to the thermodynamic limit, and exhibiting an
interesting concentration property for the distribution of values of
$\lambda$, both at equilibrium and in the non-equilibrium regime. Finally, in
Section 6.4 we assess the efficiency of the Norton approach in terms of
asymptotic variance, showing that in some cases the standard NEMD approach is
significantly outperformed.

###  6.1 Description of the numerical experiments

#### Presentation of the system.

In all experiments, we consider perturbations of the kinetic Langevin dynamics
for a Lennard–Jones fluid. The potential energy function is given by

$V_{\mathrm{LJ}}(q)=\sum_{1\leqslant i<j\leqslant
N}v\left(\|q_{i}-q_{j}\|\right),$ (83)

where $v$ is the radial function

$v(r)=4\varepsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right).$

Note that $v^{\prime}(2^{1/6}\sigma)=0$. The parameters $\varepsilon$ and
$\sigma$ modulate respectively the energy and spatial range of the
interaction. It is convenient to state results in the system of reduced units
in which $\varepsilon$, $2^{1/6}\sigma$, the mass $m$ of each particle and the
Boltzmann constant $k_{\mathrm{b}}$ are set to $1$. In fact, we consider a
slightly modified version of the potential (83), obtained by truncating the
range of $v$ as

$v_{r_{\mathrm{c}}}(r)=\left[v(r)-v(r_{\mathrm{c}})-v^{\prime}(r_{\mathrm{c}})(r-r_{\mathrm{c}})\right]\mathbbm{1}_{r\leqslant
r_{\mathrm{c}}},$

where $r_{\mathrm{c}}$ is a cutoff radius which prescribes the maximum range
of interactions, the added affine term ensuring that $v_{r_{\mathrm{c}}}$ is
$C^{1}$. We set $r_{\mathrm{c}}=2.5$ in our simulations, which were performed
using the Julia package Molly [19]. Both the Norton and NEMD simulations were
carried out using a $\mathrm{BAOAB}$ numerical splitting scheme.

#### Discrete estimators of transport coefficients.

We compute approximations of the mobility by plotting the average response as
a function of the forcing magnitude $\eta$, and fitting the slope of the
initial linear regime to obtain the transport coefficient. More precisely,
continous time estimators for the transport coefficient (7) and its Norton
analog (25) can be obtained by defining appropriate discretizations thereof.
These discretizations are given by the following estimators

$\widehat{\alpha}^{\Delta
t}_{N_{\mathrm{iter}},\eta}=\frac{1}{N_{\mathrm{iter}}\eta}\sum_{k=0}^{N_{\mathrm{iter}}-1}R(q^{k},p^{k}),\qquad\widehat{\alpha}^{\Delta
t,*}_{N_{\mathrm{iter}},r}=rN_{\mathrm{iter}}\left(\sum_{k=0}^{N_{\mathrm{iter}}-1}\lambda^{k}\right)^{-1},$
(84)

respectively in the NEMD case and the Norton case, where $\lambda^{k}$ is
defined in (82) for the Norton dynamics preserving some response $R$. More
reliable estimators can further be obtained by fitting the response profile
for several small values of $\eta$ and $r$.

In the case of mobility computations, the observable $R(q,p)=F\cdot M^{-1}p$
is used, yielding the following NEMD estimator for the mobility:

$\widehat{\alpha}^{\Delta
t}_{N_{\mathrm{iter}},\eta}=\frac{1}{N_{\mathrm{iter}}\eta}\sum_{k=0}^{N_{\mathrm{iter}}-1}F\cdot
M^{-1}p^{k}.$ (85)

For shear viscosity, discretizations of the NEMD estimator (51) for the
Fourier coefficient $U_{1}$ are used, yielding

$\widehat{U}_{1,N_{\mathrm{iter}},\eta}^{\Delta
t}=\frac{1}{NN_{\mathrm{iter}}\eta}\sum_{k=0}^{N_{\mathrm{iter}}}\sum_{j=1}^{N}\left(M^{-1}p^{k}\right)_{j,x}\exp\left(\frac{2\mathrm{i}\pi
q^{k}_{j,y}}{L_{y}}\right).$ (86)

#### Verification of the controllability condition.

We conclude this section by discussing the controllability condition (54) in
the case of mobility and shear viscosity computations.

When computing the mobility using a constant forcing, note that, using the
expression for the non-equilibrium response in (38), the controllability
condition writes

$(G\cdot F)(q)=F^{\intercal}M^{-1}F\neq 0,$

since $M^{-1}$ is a positive definite matrix, so that the controllability
condition is trivially satisfied everywhere. Moreover, the existence of strong
solutions to (52) can be straightforwardly deduced from the existence of
strong solutions to (35), since in this case
$\overline{P}_{F,G}=\overline{P}_{F,M^{-1}F}$ is a constant matrix, which
therefore preserves locally Lipschitz maps which grow at most linearly at
infinity, see for instance the discussion in [36, Section 3.3].

We examine the controllability condition in the case of a transverse shear
forcing profile in the case $M=m\mathrm{Id}$ in order to alleviate the
notation (the extension to more general mass matrices being straightforward).
Using the expression for the response (50), the controllability condition
writes

$(G\cdot F)(q)=\frac{1}{mN}\sum_{j=1}^{N}\exp\left({\frac{2\mathrm{i}\pi
q_{j,y}}{L_{y}}}\right)f(q_{j,y})\neq 0.$

The quantity $G\cdot F$ could in principle vanish. As the number of particles
is increased, this should however be rather unlikely when the first Fourier
coefficient of the forcing $f$ is non-zero. Indeed, the marginal distribution
of a single particle is uniform over the simulation cell at equilibrium (by
translation invariance). Therefore, one expects, for a large number of
particles and not too strong a forcing, that $G\cdot F$ is close to

$\frac{1}{mL_{y}}\int_{0}^{L_{y}}\exp\left({\frac{2\mathrm{i}\pi
y}{L_{y}}}\right)f(y)\,\mathrm{d}y,$

which is exactly the first Fourier coefficient of $f$, a non-zero quantity.

###  6.2 Equivalence of (non-)equilibrium responses

We begin by checking the consistency between the standard NEMD approach and
the Norton approach, by applying these two methods to the cases of mobility
and shear viscosity computations discussed in Section 4.2.

#### Mobility computations.

We begin by checking the validity of the Norton method on the simple example
of mobility. We used identical equilibrium conditions for NEMD and Norton
computations, namely a three-dimensional system of $N=1000$ particles, in a
cubic periodic domain of side length $L$, such that $\rho=N/L^{3}=0.6$, with a
temperature $T=1.25$ and a friction coefficient $\gamma=1.0$. The time step
was set to ${\Delta t}=10^{-3}$ (this choice ensures that the relative
variations of the energy of the system are of the order of 0.6 for Hamiltonian
dynamics), and simulations were performed for a minimal physical time of
$T_{\mathrm{sim}}=250,000$ in the linear response regime, and
$T_{\mathrm{sim}}=8,000$ far in the non-linear response regime. Error bars
have been omitted for the sake of readability, as they are in all cases
smaller than the size of the markers.

In Figure 2, we plot the response as a function of the forcing magnitude for
both the Norton and NEMD dynamics. Thus, the fixed quantity is plotted in the
ordinates for the Norton system, and in the abscissas for the NEMD system.
Least-squares linear regression lines are plotted, in a dotted red line for
the Norton system, and a dashed green line for the NEMD system. The slopes of
these lines, indicated in the legend, give the estimated transport
coefficient. The left plot corresponds to a color drift perturbation (40),
while the right plot corresponds to a single drift perturbation (39). We
observe that, while the agreement is almost perfect in the case of the color-
drift forcing, there is a significant discrepancy in the linear responses in
the case of a single-drift forcing. In fact, using the relation (43), it is
readily checked that the Norton estimation of the diffusion coefficient is
incorrect. We believe that for the Norton method to be valid, the microscopic
perturbation should as a general rule act on the bulk of the system. This
condition is not satisfied in the case of a single-drift forcing, which acts
on a single particle.

Figure 2: Response profiles in the linear regime for the NEMD and Norton
mobility dynamics. The slope of the linear regression lines correspond to the
estimated transport coefficient. Left: color drift. Right: single drift.

In Figure 3, we again look at the color-drift perturbed dynamics, but this
time take the system far from equilibrium, so that the response is non-linear.
Remarkably, the non-linear responses profiles for the Norton and NEMD systems
still coincide for extreme values of the forcing magnitude. This was already
formally proven in [11] in a deterministic setting.

Figure 3: Non-linear response profiles for the NEMD and Norton mobility
dynamics with color drift. The responses coincide even in the large
perturbation regime.

#### Shear viscosity computations

In Figures 4 to 6, we perform an experiment similar to the one done to
estimate the mobility, for a system subject to a longitudinal forcing
modulated in intensity by a transverse profile as described in Section 4.2,
for the three forcing profiles considered in [24]. The systems were simulated
under the same conditions as for the mobility, except for the temperature and
particle density which were respectively chosen to be $T=0.8$ and
$\rho=N/(L_{x}L_{y}L_{z})=0.7$.

More precisely, Figure 4 presents the response obtained for a piecewise linear
forcing profile $f(y)=4\mathbbm{1}_{0\leqslant
y<L_{y}/2}(y-L_{y}/4)/L_{y}+4\mathbbm{1}_{L_{y}/2\leqslant
y<L_{y}}(3L_{y}/4-y)/L_{y})$, Figure 5 for a piecewise constant forcing
profile $f(y)=1-2\mathbbm{1}_{0<y\leqslant L_{y}/2}$, and Figure 6 for a
sinusoidal forcing profile $f(y)=\sin(2y\pi/L_{y})$. In each case, the linear
response regime is plotted on the left-hand side, and the non-linear regime is
plotted on the right. Linear regression lines whose slopes are given in the
legend are superimposed to the data points. Each time, both the linear and
non-linear responses for usual NEMD and Norton dynamics perfectly agree.

#### Consistency in the thermodynamic limit.

We now turn to investigating the behavior of estimators of the shear viscosity
in the thermodynamic limit $N\to\infty$. In order for the thermodynamic limit
of the viscosity to be well-defined, the quantity $F_{1}/U_{1}-\gamma_{x}$
must scale asymptotically as $N^{-2/3}$. Indeed, as the system size increases,
the length $L_{y}$ scales as $(N/\rho)^{1/3}$, which, plugged in the
expression

$\nu=\rho\left(\frac{F_{1}}{U_{1}}-\gamma_{x}\right)\left(\frac{L_{y}}{2\pi}\right)^{2}$

implies the claimed scaling for $F_{1}/U_{1}-\gamma_{x}$. This discussion
determines how to fit the asymptotic behavior of $F_{1}/U_{1}$ for increasing
system sizes (see Figure 7).

In Figure 7, we examine the behavior of the response for increasing sizes of
the system. Estimations of $U_{1}$ were performed at roughly equivalent state
points $\eta=0.3$ and $r=0.1$ in the linear response regime using estimators
of the form (84). We find in particular that the scaling law
$F_{1}/U_{1}-\gamma_{x}\sim N^{-2/3}$ is obeyed, and that the estimates of the
prefactor $C$ coincide for the NEMD and Norton methods. The agreement between
the two methods is excellent for all values of $N$.

Figure 4: Non-equilibrium normalized Fourier response for shear viscosity
computations with a piecewise linear forcing profile.

Figure 5: Non-equilibrium normalized Fourier response for shear viscosity
computations with a piecewise constant forcing profile.

Figure 6: Non-equilibrium normalized Fourier response for shear viscosity
computations with a sinusoidal forcing profile.

###  6.3 Concentration properties in the thermodynamic limit

We next investigate the concentration rates of $R$ and $\lambda$ as a function
of the system size $N$, respectively in the canonical equilibrium and Norton
ensembles, in the case of a sinusoidal shear forcing profile. In Figure 8, we
show equilibrium distributions of $R$ and $\lambda$ for increasing values of
$N$. By ‘̀equilibrium,̈ it is meant here that the NEMD dynamics is run for
$\eta=0$ and the Norton dynamics for $r=0$. Note first that both distributions
are centered around 0 and symmetric. Moreover, they concentrate around 0 as
the system size $N$ increases. The concentration rate however differs for the
two dynamics. In the canonical ensemble, the usual concentration rate
$N^{-1}$, consistent with some form of spatial averaging of an intensive
quantity, is observed. Indeed, the rate $N^{-1}$ is the one dictated by a
central limit theorem, which makes sense on an intuitive level for a system
displaying only short-range interactions, using a decomposition of the full
system into weakly correlated subsystems, each of them contributing additively
to the total flux.

On the other hand, the concentration rate as a function of the system size in
the Norton ensemble appears to be faster. This is confirmed in Figure 9, which
shows that the second moment of $\lambda$ decays as $N^{-5/3}$. The increase
in the rate is presumably due to the coupling introduced by the constant-flux
condition, which introduces global correlations in the systems, and hence some
form of long range interactions, which are known to lead to significant
changes of phenomenology in some situations.

Figure 7: Behavior of estimators for $U_{1}/F_{1}-\gamma_{x}$ in the large $N$
limit. The simulation results, indicated by markers, are fitted by
$CN^{-\alpha}$ with $C\approx 14$ and $\alpha\approx 0.66$ with a linear
regression of $\log(U_{1}/F_{1}-\gamma_{x})$ as a function of $\log N$. The
right plot corresponds to the left plot in log-log coordinates.

Figure 8: Histograms of forcing and response in equilibrium ensembles. Left:
standard equilibrium Langevin dynamics ($\eta=0$). Right: Langevin–Norton
dynamics at zero flux ($r=0$). Figure 9: Variance of $\lambda$ and $R$ in the
Norton and canonical equilibrium ensembles.

###  6.4 Asymptotic variance

We finally assess the numerical efficiency of the Norton method. NEMD
estimators of the form (84) are ergodic averages of $R$ for trajectories of
the discrete-time Markov chain defined by the numerical scheme, and the
central limit theorem for Markov chains implies that the variance is
asymptotically $\eta^{-2}N_{\mathrm{iter}}^{-1}\sigma^{2}_{{\Delta t}}(R)$,
where $\sigma^{2}_{{\Delta t},\eta}(R)$ is the asymptotic variance of $R$
under the scheme’s stationary measure $\pi_{\Delta t,\eta}$:

$\sigma^{2}_{\Delta t,\eta}(R)=\mathrm{Var}_{\pi_{\Delta
t,\eta}}(R)+2\sum_{k=1}^{\infty}\mathrm{Cov}_{\pi_{\Delta
t,\eta}}\left(R(q^{k},p^{k})R(q^{0},p^{0})\right)=\mathrm{Var}_{\pi_{\Delta
t,\eta}}(R)\Theta_{\pi_{\Delta t,\eta}}(R),$ (87)

where

$\Theta_{\pi_{\Delta
t,\eta}}(R)=\left[1+2\sum_{k=1}^{\infty}\frac{\mathrm{Cov}_{\pi_{\Delta
t,\eta}}\left(R(q^{k},p^{k})R(q^{0},p^{0})\right)}{\mathrm{Var}_{\pi_{\Delta
t,\eta}}(R)}\right]$ (88)

is the number of correlation steps $R$ for stationary initial distribution,
and $\mathrm{Var}_{\pi_{\Delta t,\eta}}(R)$ denotes the centered second moment
of $R$ under $\pi_{\Delta t,\eta}$. It can be proved for standard NEMD
dynamics [34] that

$\sigma^{2}_{\Delta t,\eta}(R)=\sigma^{2}_{\eta}(R)/\Delta t+\mathrm{O}(1),$

where $\sigma_{\eta}^{2}(R)$ is the asymptotic variance (10) of the continuous
process.

For the Norton estimator in (84), it can be shown, using the delta method and
a computation similar to (26), that the variance of trajectory averages of $R$
is asymptotically $N_{\mathrm{iter}}^{-1}r^{4}\mathbb{E}_{\pi_{\Delta
t,r}^{*}}[\lambda]^{-4}\sigma^{2*}_{{\Delta t},r}(\lambda)$, where
$\pi_{{\Delta t},r}^{*}$ denotes the invariant probability measure of the
Markov chain $(\lambda^{n})_{n\geqslant 1}$ defined in (82), and
$\sigma^{2*}_{{\Delta t},r}(\lambda)$ is the associated asymptotic variance
for $\lambda$. Again, we write the asymptotic variance as the product of the
centered second moment under the steady-state $\mathrm{Var}_{\pi_{\Delta
t,r}^{*}}(\lambda_{0})$ and the number of correlation steps
$\Theta_{\pi_{\Delta t,r}^{*}}(\lambda)$.

As the number of particles $N$ grows to infinity, and for equivalent state
points $r=\eta\alpha$, say in the linear response regime, the numerical
results of Section 6.3 show that the contribution of the stationary centered
second moment will be asymptotically smaller for Norton estimators than for
their NEMD counterparts, owing to the fast concentration rate of $\lambda$ in
the Norton ensemble. This is reason to suspect that Norton estimators of the
transport coefficient may have lower asymptotic variance for large enough $N$,
as long as the scaling of correlation times for $\lambda$ in the Norton
ensemble does not cancel this effect. To assess whether this is the case, we
plot in Figure 10 autocorrelation functions for equilibrium trajectories of
$\lambda$ in the Norton ensemble and $R$ in the canonical ensemble, for
several values of the system size $N$, again in the case of a sinusoidal shear
forcing profile. These correspond to the summand in (88), plotted as a
function of the physical simulation time. We observe that the correlation time
is roughly independent of $N$. In fact, the autocorrelation profile itself
barely changes with increasing system sizes.

Furthermore, we observe that the correlation time is much smaller in the
Norton ensemble than in the NEMD ensemble. This suggests that even for
moderate system sizes, the asymptotic variance for estimators of the linear
response should be smaller for Norton systems than for their NEMD
counterparts, owing to smaller correlation times. We verify this intuition in
Figure 11 for a fixed system of $N=1000$ particles, in the context of shear
viscosity computations. The asymptotic variance for estimators of the Fourier
linear response $U_{1}$ is plotted as a function of the forcing magnitude, for
the various forcing profiles considered in Section 6.2. Note that the
asymptotic variances indeed scale as $\eta^{-2}$ for NEMD simulations and
$r^{-2}$ for Norton dynamics, at least for small values of $|r|$ and $|\eta|$.
Interestingly, we find that Norton estimators have lower variance in each
situation. We expect that this difference will be more pronounced for larger
systems.

Figure 10: Autocorrelation functions at equilibrium for the shear-viscosity
Fourier response $R$ (given by (50)). Left: autocorrelation function of
$\lambda$ for the Norton dynamics on $R^{-1}\\{0\\}$. Right: autocorrelation
function of $R$ in standard equilibrium simulations. Figure 11: Asymptotic
variance for estimators of the Fourier response $U_{1}$ as a function of the
perturbation magnitude, in log-log coordinates. The expected scaling line for
small values of the forcing or response is plotted in red.

##  7 Perspectives for future work

We list in this final section some of the many open questions raised by the
stochastic Norton approach, both on the theoretical and numerical sides.
Certain points are currently being tackled, but a lot remains to be done, and
we hope that this work will encourage other researchers to study Norton
dynamics. We classify issues to address in three families.

#### Analysis of the continuous dynamics.

A first task is to show that Norton dynamics are well-posed and well-behaved.
As such results will undoubtedly constitute necessary tools to investigate the
properties of Norton ensembles from a rigorous standpoint, they are of the
utmost importance. We identify the following questions:

  * $\circ$

Obtain sufficient conditions on the reference dynamics, $F$ and $R$, for the
Norton dynamics to be well-posed on arbitrarily large time horizons.

  * $\circ$

Prove the existence and uniqueness of a Norton steady-state.

  * $\circ$

Show pathwise ergodicity of the Norton dynamics.

  * $\circ$

Obtain bounds on the rate of convergence of arbitrary initial conditions to
the Norton steady-state.

In view of conditions such as (22), it is clear that not all pairs $(F,R)$
give rise to well-posed Norton dynamics. From a macroscopic point of view, the
averages of $F$ and $R$ correspond to conjugate variables, hence we expect any
reasonable microscopic realization to appropriately relate the associated
observables.

#### Theoretical questions.

We list here various physically motivated questions, which would in particular
establish that the Norton method provides a viable way to compute transport
coefficients.

  * $\circ$

Derive equivalence of ensemble results between the Norton and canonical
ensemble at equilibrium (i.e. $r=0$ and $\eta=0$, respectively), in the
thermodynamic limit.

  * $\circ$

Prove the existence of the linear response (24) for all systems of finite
size, and in fact develop a theory of linear response theory for Norton
ensembles. This would include relating the transport coefficient to
fluctuations in the Norton equilibrium ensemble, analogously to standard
Green–Kubo formulæ in the canonical ensemble. This approach was formally
carried out in [15], for the special case of electrical conductivity or
mobility, in a deterministic setting. Carrying out this development would in
particular require understanding, for small $|r|$, the relation of the Norton
equilibrium steady-state to the invariant measure of the reference process
(2), as well as to the nature of the forcing $F$ and the response $R$.

  * $\circ$

Show the consistency of the Norton approach, by proving that the transport
coefficients obtained with NEMD and Norton dynamics have the same
thermodynamic limit, for bulk forcings. This can be considered as some form of
equivalence of ensembles in a perturbative nonequilibrium setting.

  * $\circ$

Derive a result on the equivalence of ensembles for the Norton and NEMD
approaches beyond the linear response regime, for equivalent state points
(i.e. the average forcing experienced in Norton dynamics coincides with the
fixed forcing in NEMD simulations).

#### Numerical analysis of Norton dynamics.

A last set of issues to address is related to properties of Norton dynamics
which are of practical interest for the efficient and accurate computation of
transport coefficient in realistic systems. A crucial point is to derive clear
conditions under which one should prefer the Norton approach to the usual NEMD
method. Another point concerns the study of the mathematical properties of the
discretized dynamics.

  * $\circ$

Carry out a rigorous analysis of the concentration rate of $\lambda$ in the
thermodynamic limit.

  * $\circ$

Explain the shorter correlation times for Norton dynamics compared to the
corresponding NEMD dynamics. This and the previous question would ideally
yield rigorous asymptotic formulas for both the variance and correlation time
in terms of $\eta$, $r$ and $N$.

  * $\circ$

Perform the numerical analysis of splitting schemes for Norton dynamics. This
analysis is a necessary step to derive error bounds relating quantitative
properties of the continuous dynamics with their discrete analogs. Useful
results would include error estimates at the level of discrete trajectories,
or at the level of the invariant measures of the associated Markov chains.

##  Appendix A Proofs of some technical results

We gather in this appendix the proofs of various technical statements, namely
Equation (29) in Section A.1, Equation (34) in Section A.2, and Equation (43)
in Section A.3.

###  A.1 Derivation of Norton dynamics in the case of multiple constraints

We write here the proof of (29). We proceed as in Section 3.2, assuming that
the forcing process $\Lambda^{\mathbf{r}}$ may be written as an Itô process
before checking a posteriori that this assumption is justified. More
precisely, we write

$\mathrm{d}\Lambda_{t}^{\mathbf{r}}=\mathrm{d}\widetilde{\Lambda}_{t}^{\mathbf{r}}+\lambda^{\mathbf{r}}_{t}\,\mathrm{d}t,$

where $\widetilde{\Lambda}^{\mathbf{r}}$ is the martingale part of the Itô
decomposition of $\Lambda^{\mathbf{r}}$. Applying Itô’s formula to the
constant-response constraint
$R(Y^{\mathbf{r}}_{t})=R(Y^{\mathbf{r}}_{0})=\mathbf{r}$ yields

$\nabla
R(Y^{\mathbf{r}}_{t})^{\intercal}\mathrm{d}Y^{\mathbf{r}}_{t}+\frac{1}{2}\nabla^{2}R(Y^{\mathbf{r}}_{t}):\mathrm{d}\langle
Y^{\mathbf{r}}\rangle_{t}=0.$ (89)

Identifying martingale increments, the uniqueness of the Itô decomposition
implies

$\mathrm{d}\widetilde{\Lambda}_{t}^{\mathbf{r}}=-\left[\nabla
R(Y_{t}^{\mathbf{r}})^{\intercal}F(Y_{t}^{\mathbf{r}})\right]^{-1}\nabla
R(Y_{t}^{\mathbf{r}})^{\intercal}\sigma(Y_{t}^{\mathbf{r}})\,\mathrm{d}W_{t},$

which in turns allows to compute the quadratic covariation term

$\mathrm{d}\langle Y^{\mathbf{r}}\rangle_{t}=\left[\overline{P}_{F,\nabla
R}\sigma\sigma^{\intercal}\overline{P}_{F,\nabla
R}^{\intercal}\right](Y_{t}^{\mathbf{r}})\,\mathrm{d}t=\Pi_{F,\nabla
R,\sigma}(Y_{t}^{\mathbf{r}})\,\mathrm{d}t.$

Using this expression in (89) and the uniqueness of the Itô decomposition once
again to identify bounded-variation parts, we recover

$\lambda_{t}^{\mathbf{r}}=-(\nabla R^{\intercal}F)^{-1}\left[\nabla
R^{\intercal}b+\frac{1}{2}\left(\nabla^{2}R^{\intercal}:\Pi_{F,\nabla
R,\sigma}\right)\right](Y_{t}^{\mathbf{r}}),$

which yields (29) upon substituting the expression for
$\mathrm{d}\Lambda_{t}^{\mathbf{r}}$.

###  A.2 Derivation of the Norton dynamics in the case of time-dependent
constraints

To prove (34), we once again proceed as in Section 3.2, writing the Itô
decomposition for the forcing process as

$\mathrm{d}\Lambda_{t}^{\mathcal{R}}=\mathrm{d}\widetilde{\Lambda}_{t}^{\mathcal{R}}+\lambda^{\mathcal{R}}_{t}\,\mathrm{d}t.$

Applying Itô’s formula to the time-dependent constraint (34), we get

$\nabla
R(Y^{\mathcal{R}}_{t})^{\intercal}\mathrm{d}Y^{\mathcal{R}}_{t}+\frac{1}{2}\nabla^{2}R(Y^{\mathcal{R}}_{t}):\mathrm{d}\langle
Y^{\mathcal{R}}\rangle_{t}=\overline{r}_{t}\,\mathrm{d}t+\widetilde{r}_{t}\,\mathrm{d}B_{t},$

which, by uniqueness of the Itô decomposition, leads to

$\mathrm{d}\widetilde{\Lambda}_{t}^{\mathcal{R}}=\frac{\widetilde{r}_{t}\,\mathrm{d}B_{t}-\nabla
R(Y_{t}^{\mathcal{R}})\cdot\sigma(Y_{t}^{\mathcal{R}})\,\mathrm{d}W_{t}}{\nabla
R(Y_{t}^{\mathcal{R}})\cdot F(Y_{t}^{\mathcal{R}})},$

allowing to compute the covariation bracket as

$\mathrm{d}\left\langle Y^{\mathcal{R}}\right\rangle_{t}=\Pi_{F,\nabla
R,\sigma}(Y^{\mathcal{R}}_{t})\,\mathrm{d}t+\frac{F(Y^{\mathcal{R}}_{t})\otimes
F(Y^{\mathcal{R}}_{t})}{\left(\nabla R(Y^{\mathcal{R}}_{t})\cdot
F(Y^{\mathcal{R}}_{t})\right)^{2}}\widetilde{r}_{t}^{2}\,\mathrm{d}t.$

One can then proceed to identify bounded-variation increments, and isolate
$\lambda^{\mathcal{R}}_{t}$:

$\lambda^{\mathcal{R}}_{t}=\frac{\displaystyle{\overline{r}_{t}-\frac{1}{2}\nabla^{2}R(Y^{\mathcal{R}}_{t}):\left[\Pi_{F,\nabla
R,\sigma}(Y^{\mathcal{R}}_{t})+\frac{F(Y^{\mathcal{R}}_{t})\otimes
F(Y^{\mathcal{R}}_{t})}{\left(\nabla R(Y^{\mathcal{R}}_{t})\cdot
F(Y^{\mathcal{R}}_{t})\right)^{2}}\widetilde{r}_{t}^{2}\right]-\nabla
R(Y^{\mathcal{R}}_{t})\cdot b(Y^{\mathcal{R}}_{t})}}{\nabla
R(Y^{\mathcal{R}}_{t})\cdot F(Y^{\mathcal{R}}_{t})}.$

Substituting in the expression for the forcing process finally yields (34).

###  A.3 Derivation of the relation between color/single drift linear
responses

The proof of (43) is taken from unpublished notes by Julien Roussel, see the
PhD thesis [39]. We assume for the ease of presentation that $M=m\mathrm{Id}$
and $d=3$, and denote by $\mu$ the equilibrium measure (36). The dynamics
under consideration is the standard Langevin dynamics, namely (35) with
$\eta=0$. Let us define, for $1\leqslant i,j\leqslant N$,

$c_{ij}=\frac{\beta}{m^{2}}\int_{0}^{\infty}\mathbb{E}_{\mu}[p_{i,x,t}p_{j,x,0}]\,\mathrm{d}t.$

In view of (42), upon summing over $i$ the longitudinal $p$-components of the
SDE (35) and integrating in time,

$\displaystyle\sum_{i=1}^{N}p_{i,x,t}$
$\displaystyle=\sum_{i=1}^{N}\left[p_{i,x,0}+\int_{0}^{t}\left(-\frac{\partial}{\partial
q_{i,x}}V(q_{s})-\frac{\gamma}{m}p_{i,x,s}\,\mathrm{d}s+\sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_{i,x,s}\right)\right]$
$\displaystyle=\sum_{i=1}^{N}\left[p_{i,x,0}-\frac{\gamma}{m}\int_{0}^{t}p_{i,x,s}\,\mathrm{d}s+\sqrt{\frac{2\gamma}{\beta}}W_{i,x,t}\right].$

Multiplying by $p_{1,x,0}$ and taking the expectation over trajectories
started from canonical initial conditions (so that the Brownian terms
vanishes), it follows that

$\mathbb{E}_{\mu}\left[\left(\sum_{i=1}^{N}p_{i,x,t}\right)p_{1,x,0}\right]=\sum_{i=1}^{N}\left(\mathbb{E}_{\mu}\left[p_{i,x,0}p_{1,x,0}\right]-\frac{\gamma}{m}\int_{0}^{t}\mathbb{E}_{\mu}\left[p_{i,x,s}p_{1,x,0}\right]\,\mathrm{d}s\right).$
(90)

By the decay properties of the evolution semigroup (obtained for instance by
hypocoercive approaches, see e.g. [47], as well as the introduction of [4] for
a review), the left-hand side converges to $0$ as $t\to\infty$, while the
integral is well-defined. Since $p_{0}$ has diagonal covariance with respect
to $\mu$, we get

$\sum_{i=1}^{N}\frac{\gamma}{m}\int_{0}^{\infty}\mathbb{E}_{\mu}\left[p_{i,x,s}p_{1,x,0}\right]\,\mathrm{d}s=\mathbb{E}_{\mu}[p_{1,x,0}^{2}]=\frac{m}{\beta}.$
(91)

Equivalently,

$\sum_{i=1}^{N}c_{i1}=\frac{1}{\gamma}.$ (92)

Using the indistinguishability of the particles,

$\forall\,1\leqslant i,j\leqslant N,\qquad c_{ii}=c_{11},\qquad
c_{ij}=c_{12}.$ (93)

We can therefore rewrite (92) as

$c_{12}=\frac{1}{N-1}\left(\frac{1}{\gamma}-c_{11}\right).$ (94)

This equality can be used to relate the linear responses of the single drift
and the color drift. By the Green–Kubo formula [38], the transport coefficient
for the single drift is given by

$\alpha_{F_{\mathrm{S}}}=c_{11}.$

For the color drift, we expand, using the Green–Kubo formula,

$\alpha_{F_{\mathrm{C}}}=\frac{\beta}{m^{2}}\int_{0}^{\infty}\mathbb{E}_{\mu}[\left(F_{\mathrm{C}}\cdot
p_{t}\right)\left(F_{\mathrm{C}}\cdot
p_{0}\right)]\,\mathrm{d}t=\frac{1}{N}\left(\sum_{i=1}^{N}c_{ii}+2\sum_{1\leqslant
i<j\leqslant N}(-1)^{i+j}c_{ij}\right)$

In view of (93), and using

$\sum_{1\leqslant i<j\leqslant
N}(-1)^{i+j}=-\left\lfloor\frac{N}{2}\right\rfloor,$

which is easily seen by induction, we get

$\alpha_{F_{\mathrm{C}}}=c_{11}-\frac{2\lfloor
N/2\rfloor}{N(N-1)}\left(\frac{1}{\gamma}-c_{11}\right),$ (95)

which is the claimed identity.

#### Acknowledgments.

The authors thank Gary Morris, Julien Roussel and Giovanni Ciccotti for
stimulating discussions, and the members of the project ANR SINEQ (in
particular Benjamin Jourdain, Tony Lelièvre, Stefano Olla and Mathias Rousset)
for valuable feedback on preliminary results presented at the annual meeting
of the project. This work received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (project EMC2, grant agreement No 810367), and from the Agence
Nationale de la Recherche, under grants ANR-19-CE40-0010-01 (QuAMProcs) and
ANR-21-CE40-0006 (SINEQ).

## References

  * [1] M.P. Allen and D.J. Tildesley.  Computer Simulation of Liquids.  Oxford University Press, 2017.
  * [2] H. Andersen.  RATTLE: A “velocity” version of the SHAKE algorithm for molecular dynamics calculations.  Journal of Computational Physics, 52(1):24–34, 1983.
  * [3] R. Balian.  From Microphysics to Macrophysics.  Springer, 1991.
  * [4] E. Bernard, M. Fathi, A. Levitt, and G. Stoltz.  Hypocoercivity with Schur complements.  Annales Henri Lebesgue, 5:523–557, 2022.
  * [5] M. Christian and N. Karel.  Rigorous meaning of McLennan ensembles.  Journal of Mathematical Physics, 51(1):015219, 2010.
  * [6] G. Ciccotti and G. Jacucci.  Direct computation of dynamical response by molecular dynamics: The mobility of a charged Lennard-Jones particle.  Phys. Rev. Lett., 35:789–792, 1975.
  * [7] G. Ciccotti, R. Kapral, and A. Sergi.  Non-equilibrium molecular dynamics.  Handbook of Materials Modeling, 1:1–17, 2005.
  * [8] G. Ciccotti, T. Lelièvre, and E. Vanden-Eijnden.  Projection of diffusions on submanifolds: Application to mean force computation.  Comm. Pure Appl. Math., 61(3):371–408, 2008.
  * [9] S. Darshan, A. Eberle, and G. Stoltz.  Sticky coupling as a control variate for sensitivity analysis, 2023.  In preparation.
  * [10] Matthew Dobson.  Periodic boundary conditions for long-time nonequilibrium molecular dynamics simulations of incompressible flows.  The Journal of Chemical Physics, 141(18):184103, 2014.
  * [11] D.J. Evans.  The equivalence of Norton and Thévenin ensembles.  Molecular Physics, 80:221–224, 1993.
  * [12] D.J. Evans and J.F. Ely.  Viscous flow in the stress ensemble.  Molecular Physics, 59:1043–1048, 1986.
  * [13] D.J. Evans, W.G. Hoover, B.H. Failor, B. Moran, and A.J.C. Ladd.  Nonequilibrium molecular dynamics via Gauss’s principle of least constraint.  Physical Review A, 28(2):1016–1021, 1983.
  * [14] D.J. Evans and G. Morriss.  Statistical Mechanics of Nonequilibrium Liquids.  Cambridge University Press, 2008.
  * [15] D.J. Evans and G.P. Morriss.  Equilibrium-fluctuation expression for the resistance of a Norton circuit.  Physical Review A, 31(6):3817–3819, 1985.
  * [16] D. Frenkel and B. Smit.  Understanding Molecular Simulation: from Algorithms to Applications.  Elsevier, 2001.
  * [17] E.M. Gosling, I.R. McDonald, and K. Singer.  On the calculation by molecular dynamics of the shear viscosity of a simple fluid.  Molecular Physics, 26:1475–1484, 1973.
  * [18] M.S. Green.  Markoff random processes and the statistical mechanics of time‐dependent phenomena. II. Irreversible processes in fluids.  The Journal of Chemical Physics, 22(3):1281–1295, 1954.
  * [19] J. Greener.  Molly.  _Github repository: github.com/JuliaMolSim/Molly.jl_, 2020.
  * [20] E. Hairer, G. Wanner, and C. Lubich.  Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations.  Springer, 2006.
  * [21] M. Hairer and A.J. Majda.  A simple framework to justify linear response theory.  Nonlinearity, 23(4):909–922, 2010.
  * [22] W.G. Hoover, H.A. Posch, and L.W. Campbell.  Thermal heat reservoirs via Gauss’ principle of least constraint: Dissipation, chaos, and phase‐space dimensionality loss in one‐dimensional chains.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 3(3):325–332, 1993.
  * [23] L. Hörmander.  Hypoelliptic second order differential equations.  Acta Mathematica, 119:147–171, 1967.
  * [24] R. Joubaud and G. Stoltz.  Nonequilibrium shear viscosity computations with Langevin dynamics.  Multiscale Modeling & Simulation, 10:191–216, 2012.
  * [25] W. Kliemann.  Recurrence and invariant measures for degenerate diffusions.  The Annals of Probability, 15(2):690–107, 1987.
  * [26] T. Komorowski, J. Lebowitz, S. Olla, and M. Simon.  On the conversion of work into heat: microscopic models and macroscopic equations.  _arXiv preprint 2212.00093_, 2022.
  * [27] B. Leimkuhler and C. Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Applied Mathematics Research eXpress, 2013(1):34–56, 2013.
  * [28] B. Leimkuhler and C. Matthews.  Robust and efficient configurational molecular sampling via Langevin dynamics.  The Journal of Chemical Physics, 138(17):174102, 2013.
  * [29] B. Leimkuhler and C. Matthews.  Molecular Dynamics.  Springer, 2015.
  * [30] B. Leimkuhler, C. Matthews, and G. Stoltz.  The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics.  IMA J. Numer. Anal., 36(1):13–79, 2016.
  * [31] B. Leimkuhler and S. Reich.  Simulating Hamiltonian Dynamics.  Cambridge University Press, 2004.
  * [32] T. Lelièvre, M. Rousset, and G. Stoltz.  Free Energy Computations - A Mathematical Perspective.  Imperial College Press, 2010.
  * [33] T. Lelièvre, M. Rousset, and G. Stoltz.  Langevin dynamics with constraints and computation of free energy differences.  Mathematics of Computation, 81:2071–2125, 2012.
  * [34] T. Lelièvre and G. Stoltz.  Partial differential equations and stochastic methods in molecular dynamics.  Acta Numerica, 25:681–880, 2016.
  * [35] N. Metropolis, A. W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics, 21(6):1087–1092, 1953.
  * [36] G. A. Pavliotis.  Stochastic Processes and Applications, volume 60 of Texts in Applied Mathematics.  Springer, New York, 2014.
  * [37] L. Rey-Bellet.  Ergodic properties of Markov processes.  In Open Quantum Systems II: The Markovian Approach, Lecture Notes in Mathematics, pages 1–39. Springer, 2006.
  * [38] H. Rodenhausen.  Einstein’s relation between diffusion constant and mobility for a diffusion model.  J. Statist. Phys., 55(5-6):1065–1088, 1989.
  * [39] J. Roussel.  Theoretical and Numerical Analysis of Non-Reversible Dynamics in Computational Statistical Physics.  PhD thesis, Université Paris-Est, 2018.
  * [40] J.P. Ryckaert, G. Ciccotti, and H. Berendsen.  Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.  Journal of Computational Physics, 23(3):327–341, 1977.
  * [41] R. Spacek and G. Stoltz.  Extending the regime of linear response with synthetic forcings.  _arXiv preprint 2303.03551_, 2023.
  * [42] G. Stoltz.  Error estimates and variance reduction for nonequilibrium stochastic dynamics, 2022.  _arXiv preprint 2211.10717_, 2022.
  * [43] B. D. Todd and P. J. Daivis.  Homogeneous non-quilibrium molecular dynamics simulations of viscous flow: Techniques and applications.  Mol. Sim., 33(3):189–229, 2007.
  * [44] B. D. Todd and P. J. Daivis.  Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications.  Cambridge University Press, 2017.
  * [45] M. Tuckerman.  Statistical Mechanics: Theory and Molecular Simulation.  Oxford University Press, 2010.
  * [46] A. Van der Vaart.  Asymptotic Statistics, volume 3.  Cambridge University Press, 2000.
  * [47] C. Villani.  Hypocoercivity.  Mem. Amer. Math. Soc., 202(950), 2009.